.  ~  it.: 


J  ■■DmTffiBWQti"  MWS -A 


prevfed-  lor  public  release; 
Distribution  Unlimited; 


BEST 

AVAILABLE  COPY 


UNCLASSIFIED 


Security  Classification 


3200.8  (Att  1  to  Enel  l) 
Mar  7,  66 


DOCUMENT  CONTROL  DATA  .  R  &  D 


i 


1.  ORICINATINC  ACTIVITY  (ClipMM  MllMIJ 

Woodward-Lundgren  §  Associates 

Unclassified 

IS.  CROUP 

J.  REPORT  TITLE 

A  Theoretical  Method  for  Evaluating  Stability  of  Openings  in  Rock 

«.  DESCRIPTIVE  MOVE*  (7>p«  ./  fpert  an.  lnclu.fr.  Mu) 

Final  Report  -  March  12.  1971  -  Marrh  1?  icm 

Chin-Yung  Chang  and  Keshavan  Nair 

April  12,  1972 

7 A.  TOTAL  NO.  OF  PACES  1b.  NO.  OF  REFS 

144  31 

M.. CONTRACT  OR  BRANT  NO. 

H0219046 

b.  PROJECT  NO. 

c. 

A 

SA.*  ORIGINATOR**  REPORT  NUMBER!*! 

None 

M.  OTHER  REPORT  NOIs)  (Any  ath.r  nunMr.  Iliaf  may  6.  attlgntd 

Ml  report) 

None 

10.  DISTRIBUTION  STATEMENT 

Distribution  of  this  document  is  unlimited. 

II.  VumpluiinYary  notes 

It.  SPONSORING  MILITARY  ACTIVITY 

Advanced  Research  Projects 
Agency,  Washington,  DC  20301 

N 


IS.  ABITAACT 


The  purposes  of  this  study. are  to  evaluate  the  reliability  of 
the  available  finite  element  techniques  for  the  solution  of 
plane  problems  to  predict  the  stresses,  strains,  and  deforma¬ 
tions  in  a  rock  mass  surrounding  an  excavation.  The  approach  to 
this  study  is  divided  into  two  phases  j  developing  a  general 

analytical  procedure  (i.e.,  a  finite  element  computer  program 
with  wide  capabilities)  for  determining  the  mechanical  state  in  aj 
rock  mass,  and  Analysis  of  case  histories  to  compare  pre¬ 
dicted  and  observed  values. ('This  report  discusses  essential 
features  of  finite  element  tednniques  for  modelling  rock  behavior] 
and  describes  the  development  and  verification  of  a  single  gen¬ 
eral  computer  program  which  is  capable  of  modelling  joints, 
faults,  bedding  planes  or  other  geologic  discontinuities,  "no 
tension"  properties  and  elasto-plastic  behavior  of  rock  masses. 
Using  the  computational  techniques  developed  in  Phase  1,  Morrow 
Point  Underground  Powerplant  Excavation  and  t.  laboratory  model  of 
unlined  openings  in  a  rock-like  material  were  analyzed.  It  was 
shown  tliat  the  results  of  the  analysis  provided  a  good  qualita¬ 
tive  basis  for  predicting  the  performance  of  the  model  and  the 
field  study. _  • 


DD  ££.1473 


UNCLASSIFIED 

Security  Classification 


8 


UNCLASSIFIED' 

Security  Classification 


Underground  opening 
No  tension  analysis 
Joint  perturbation  analysis 
Elasto-plastic  analysis 
Stress  transfer  technique 
Plane  strain 

Morrow  Point  Underground  Powerplant 
Model  tests  of  unlined  openings 


3200.8  (Att  1  to  Enel  l) 
Mar  7>  66 


9 


UNCLASSIFIED 

Security  Classification 


FINAL  REPORT 


MARCH  12,  1971  -  MARCH  12,  1972 


ARPA  Order  Number:  1579,  Amend.  2  Contract  Number:  H0210046 

Program  Code  Number:  1F10  Principal  Investigator: 

K.  Nair 

Telephone  Number:  (415)  444-1256 

Project  Scientist  or  Engineer: 

C-Y  Chang,  K.  Nair 

Telephone  Number:  (415)  444-1256 

Short  Title  of  Work: 

A  Theoretical  Method  for 
Evaluating  Stability  of 
Openings  in  Rock 

Amount  of  Contract: 

$34,829 


Effective  Date  of  Contract: 
March  12,  1971 

Contract  Expiration  Date: 
March  12,  1972 


Name  of  Contractor: 

Woodward-Lundgren  8  Associates 


This  research  was  supported  by  the  Advanced  Research  Projects 
Agency  of  the  Department  of  Defense  and  was  monitored  by  Bureau  of 
Mines  under  Contract  Number  H0210046. 


WOODWARD-LUNDGREN  &  ASSOCIATES 


-1- 


FINAL  REPORT  SUMMARY 


Contract  Objectives 

A  major  consideration  in  the  design  of  excavations  in  rock  is 
the  evaluation  of  the  structural  stability  of  the  opening.  An 
essential  step  in  this  evaluation  is  the  determination  of  the 
mechanical  state  in  the  rock  mass  in  the  vicinity  of  an  opening 
through  the  formulation  and  solution  of  boundary  value  problems. 
The  finite  element  method  has  proven  to  be  a  very  powerful  tool 
in  the  solution  of  boundary  value  problems.  Within  the  general 
framework  of  the  finite  element  method,  procedures  have  been 
developed  to  model  the  response  of  certain  classes  of  joints, 
cracks  and  fissures  in  the  rock  mass,  the  inability  of  rocks  to 
withstand  tension,  localized  yielding  of  rock  due  to  stress  con¬ 
centrations  and  the  time -dependent  (creep)  response  of  rock. 
While  preliminary  analysis  has  indicated  that  these  procedures 
have  the  potential  for  predicting  performance,  their  use  in 
design  will  be  limited  unless  their  reliability  is  established. 
The  objective  of  this  contract  is  to  evaluate  the  ability  of 
certain  available  finite  element  techniques  for  the  solution  of 
plane  problems  to  predict  the  stresses,  strains  and  deformations 
in  a  rock  mass  surrounding  an  excavation. 

General  Approach  and  Technical  Results 

The  approach  to  this  study  can  be  divided  into  two  phases: 

(1)  Developing  computational  procedures  which  model  the 


WOODWARD-LUNDGREN  6,  ASSOCIATES 


- 11  - 


essential  features  known  to  be  present  in  a  rock  mass,  for  deter¬ 
mining  its  mechanical  state,  pages  1  to  60,  and  (2)  Analysis  of 
case  histories  and  comparison  of  predicted  and  measured  performance, 
pages  61  to  102. 

The  essential  features  of  the  various  finite  element  techniques 
to  conduct  the  following  analyses  were  reviewed:  (I)  "No-tension" 
Analysis,  (II)  Joint  Perturbation  Analysis,  (III)  Elastic -Plastic 
Analysis,  and- (IV)  Time-dependent  (Creep)  Analysis.  Because  the 
time -dependent  (creep)  response  of  the  great  majority  of  rocks  is 
not  as  significant  as  the  time-independent  response,  and  further¬ 
more,  since  the  computational  procedures  for  time  - independent  and 
time -dependent  problems  are  significantly  different,  it  was  decided 
that  only  time -independent  analyses  for  plane  problems  would  be 
considered  in  this  study.  It  was  found  that  different  computa¬ 
tional  procedures  were  utilized  to  conduct  the  above-mentioned 
time  - independent  analyses.  In  order  to  develop  a  single  computer 
program  which  could  be  used  for  conducting  analyses  which  include 
the  features  in  (I),  (II),  and  (III),  it  was  necessary  to  formu¬ 
late  a  consistent  computational  procedure  and  develop  a  program 
on  the  basis  of  such  a  procedure.  The  "initial  stress"  (stress 
transfer)  technique  presented  by  Zienkiewicz  and  his  co-workers 
was  used  and  an  operational  finite  element  program  with  the 
capability  of  conducting  a  "no  tension,"  joint  perturbation  and 
elastic -plastic  analysis  was  developed.  Various  illustrative 
examples  were  solved  using  the  combined  program  and  the  results 
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compared  where  possible  with  results  published  by  other 
investigators. 

A  number  of  case  history  and  laboratory  studies  were  considered 
for  analysis  using  the  computational  techniques  developed  in 
Phase  I.  One  case  history  study  (Morrow  Point  Power  Plant)  and 
one  laboratory  model  were  analyzed,  and  the  results  of  the  anal¬ 
ysis  provided  a  good  qualitative  basis  for  predicting  the  per¬ 
formance  of  the  model  and  the  field  study.  Quantitative  estimates 
were  of  the  right  magnitude.  The  major  difficulty  was  in  ideal¬ 
izing  the  problem.  Since  there  are  always  features  in  the  actual 
problem  which  cannot  be  accurately  modelled,  it  is  necessary  to 
introduce  various  approximations.  In  actual  practice,  more  than 
one  idealization  of  a  problem  will  have  to  be  considered.  It  is 
believed  that  the  techniques  of  analysis  provide  the  basis  for 
identifying  problems  and  suggest  possible  design  solutions. 

POD  Implications 

The  evaluation  of  the  structural  stability  of  underground  openings, 
ground  support  structures,  and  other  facilities  is  an  essential 
step  both  in  the  design  and  in  the  survivability/vulnerability 
assessment  of  underground  structures  and  the  weapon  systems.  The 
computer  code  developed  for  this  contract  is  a  tool  that  can  be 
utilized  to  evaluate  the  stability  of  unsupported  tunnel  openings 
sited  in  rock  masses  where  the  rock  mass  behavior  is  dominated 
by  block  slippage  along  discrete  joint  planes,  or  a  global  inability 
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of  the  rock  mass  to  resist  tensile  stress,  or  elastic- 
plastic  behavior  of  the  rock  mass,  or  any  combination  of  the 
three  rock  mass  physical  characteristics.  Under  certain  cir¬ 
cumstances,  simple  ground  support  structures  may  be  modeled 
as  tunnel  liners,  although  this  feature  of  the  computer  program 
has  not  been  exploited  in  this  contract. 

Considerations  for  Further  Research 

Further  research  should  consider  increasing  the  capability  of 
the  computational  procedures  to  include  construction  sequence, 
rock  reinforcement  schemes,  and  commonly  used  steel  supports. 
Computational  techniques  to  predict  and  examine  the  influence 
of  progressive  failure  should  be  developed.  Efforts  should  be 
directed  towards  incorporating  the  work  done  by  other  investi¬ 
gators  to  extend  the  analytical  capabilities  developed  under 
this  contract. 

Additional  case  studies  should  be  investigated  to  further  eval¬ 
uate  the  capabilities  of  the  computer  program.  Research  in  the 
determination  of  in-situ  rock  mass  properties  and  conditions  is 
necessary  to  fully  utilize  the  capabilities  of  existing  analytical 
techniques . 
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The  development  of  theoretically  sound  methods  for  designing 
excavations  in  rock  is  of  significance  because  of  the  increased 
use  of  underground  facilities  for  civilian  and  military  purposes. 

A  major  consideration  in  the  design  of  excavations  in  rock  is  the 
evaluation  of  the  structural  stability  of  the  opening.  If  a 
mechanistic  approach  is  taken  for  the  evaluation  of  the  struc¬ 
tural  stability,  then  an  essential  step  in  the  evaluation  is  the 
determination- of  the  mechanical  state  (i.e.,  stresses,  strains, 
and  deformation)  in  the  rock  mass  in  the  vicinity  of  the  excava¬ 
tion.  In  recent  years,  new  technology  has  been  developed  which 
has  the  potential  of  predicting,  with  greater  accuracy,  than 
heretofore  possible,  the  mechanical  state  in  the  rock  mass. 

This  would  place  the  evaluation  of  the  structural  stability  of 
an  opening  in  rock  on  a  theoretically  sound  basis.  This  new 
technology  has  been  primarily  in  the  area  of  numerical  methods 
for  the  solution  of  boundary  value  problems.  However,  the 
application  of  these  methods  to  practical  design  problems  has 
been  very  limited.  The  primary  reason  for  this  is  the  designer's 
lack  of  confidence  in  the  ability  of  theoretical  methods  to  assist 
him  in  analyzing  practical  design  problems. 

Of  the  numerical  methods  developed,  the  finite  element  method 
has  proven  to  be  the  most  powerful  for  the  solution  of  boundary 
value  problems  (Clough,  1960,  Wilson,  1963,  1965,  Zienkiewicz, 
1967).  The  initial  application  rf  the  finite  element  method  to 
rock  mechanics  problems  used  techniques  (computer  programs)  that 
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had  been  developed  for  the  structural  analysis  of  linear 
elastic  continuous  structures.  When  the  results  obtained 
from  using  these  techniques  were  compared  with  field  results, 
they  were  found  to  be  inadequate  for  predicting  the  response 
of  rock  masses  in  the  vicinity  of  an  excavation.  Typical 
results  from  such  studies  have  been  presented  by  Judd  and 
Perloff  (1971).  It  was  recognized  that  one  of  the  major 
reasons  for  the  discrepancies  between  computed  and  observed 
behavior  was  the  inability  of  the  computational  techniques 
utilized  to  include  the  natural  geologic  discontinuities  that 
were  usually  prevalent  in  a  rock  mass.  These  discontinuities 
could  exist  in  the  rock  formation  prior  to  the  excavation  or 
could  be  a  result  of  the  induced  stresses  due  to  the  excavation 
causing  failure  in  the  rock.  Methods  of  analysis  which  have 
the  capability  of  modelling  joints  and  other  forms  of  discon¬ 
tinuities  that  commonly  exist  in  rock  masses  have  been  developed 
by  Goodman,  Taylor,  and  Brekke  (1968),  and  Zienkiewicz,  et  al . 
(1970).  Methods  for  stress  analysis  in  a  rock  mass  which  can¬ 
not  sustain  tensions  due  to  the  presence  of  cracks  and  fissures 
have  also  been  developed,  Zienkiewicz,  et  al .  (1968).  In  addi¬ 
tion,  methods  for  analyzing  localized  failure  in  a  rock  mass 
due  to  yielding  have  been  presented  by  Reyes  and  Deere  (1966)  . 
Approximate  techniques  for  incorporating  nonlinear  time-dependent 
material  properties  have  been  developed  by  Nair  and  Boresi  (1970). 
Preliminary  analysis  has  indicated  that  these  new  methods  of 
analysis  using  finite  element  techniques  have  the  potential  of 
predicting  performance  with  improved  accuracy.  However,  the 
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istic  e.g.,  no  tension  or  joints,  are  not  adequate  to  study  the 
general  case  where  all  these  factors  may  be  present.  Furthermore, 
there  has  been  very  limited  verification  of  these  techniques  on 
the  basis  of  comparison  with  measured  field  performance.  Without 
field  verification,  the  use  of  these  new  techniques  in  the  design 
of  excavations  in  rock  will  remain  limited.  Therefore,  to  develop 
a  theoretically  sound  method  for  designing  excavations  in  rock, 
which  will  be' used  in  practice,  an  essential  step  is  to  establish 
the  reliability  of  the  available  methods  of  stress  analysis  in 
predicting  the  behavior  of  rock  masses. 

PURPOSE 

The  purpose  of  this  study  is  to  evaluate  the  ability  of  the 
available  finite  element  techniques  for  the  solution  of  plane 
problems  to  predict  the  stresses,  strains,  and  deformations  in 
a  rock  mass  surrounding  an  excavation. 

GENERAL  APPROACH 

The  approach  to  this  study  can  be  divided  into  two  phases: 

(1)  Development  of  a  general  analytical  procedure  (i.e.,  a 
finite  element  computer  program  with  wide  capabilities)  for 
determining  the  mechanical  state  in  a  rock  mass,  and  (2)  Anal¬ 
ysis  of  case  histories  to  compare  predicted  and  observed  values. 

Phase  I  consisted  of  three  steps: 

a.  Review  of  pertinent  available  finite  element  programs 
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general  program. 


c.  Use  of  examples  to  illustrate  the  use  of  the  program. 

Phase  2  consisted  of  analyzing  selected  case  histories  and  com¬ 
paring  observed  and  predicted  behavior. 

ESSENTIAL  FEATURES  OF  FINITE  ELEMENT  TECHNIQUES  FOR  MODELLING 
ROCK  BEHAVIOR  (Phase  la) 

Various  techniques  have  been  utilized  in  conjunction  with  the 
finite  element  method  to  include  nonlinear  and  time-dependent 
properties,  elasto -plastic  yielding,  and  geologic  discontinuities 
in  the  stress  analysis  of  rock  masses.  These  techniques  were 
evaluated  for  the  purpose  of  developing  a  general  computer  pro¬ 
gram  for  plane  problems  which  can  incorporate  the  significant 
features  that  exist  in  rock  masses.  The  available  techniques 
use  different  computational  techniques  within  the  general  frame¬ 
work  of  the  finite  element  method  to  develop  solutions  to  bound¬ 
ary  value  problems.  The  basic  concepts  of  the  finite  element 
method  have  been  discussed  extensively  in  the  literature,  e.g., 
Clough  (1960),  Wilson  (1963,  1964)  and  Zienkiewicz  (1967).  These 
will  not  be  repeated  here.  Those  aspects  which  are  considered 
special  modifications  for  rock  mechanics  problems  are  discussed 
in  general  terms  with  respect  to  the  computational  techniques 
in  the  subsequent  paragraphs.  These  are  considered  in  four  cate¬ 
gories  (1)  "No  Tension"  analysis,  (2)  Joint  Perturbation  Analysis, 
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(3)  Elasto-Plastic  Analysis,  and.  (4)  lime  Dependent  Analysis. 
Detailed  explanation  of  these  techniques  is  provided  in  a 
subsequent  section  where  the  basis  of  the  developed  computer 
program  is  described. 

1 .  "No  Tension"  Analysis 

When  numerous  cracks  and  fissures  are  present  in  a  rock  mass, 
it  has  been  assumed  that  the  rock  is  incapable  of  withstanding 
tensile  stresses.  A  procedure  for  modelling  this  nonlinear 
behavior  has  been  presented  by  Zienkiewicz,  et  al .  (1968).  This 
method  is  called  "no  tension"  or  "strqss  transfer"  analysis. 

This  method  is  composed  of  the  following  four  essential  steps: 

(a)  A  linear  elastic  solution  to  the  problem  is  first 
obtained.  The  induced  changes  in  stress  are  added 
to  the  initial  stresses  and  the  principal  stresses 
computed . 

(b)  Those  elements  in  which  tensile  stresses  are  present 
are  identified.  As  the  material  is  assumed  incapable 
of  sustaining  them,  the  calculated  tensile  principal 
stresses  are  eliminated  without  permitting  any  point 
in  the  structure  to  displace.  In  order  to  maintain 
equilibrium,  equivalent  (balancing)  nodal  forces  are 
calculated  and  temporarily  applied  to  the  structure. 
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(c)  The  elastic  analysis  is  repeated  U-e-> 

rium  equations  are  satisfied)  to  remove  the  balancing 
nodal  forces  and  the  check  for  tensile  stresses  is 
repeated . 

(d)  If  at  the  end  of  step  (c)  tensile  stresses  are  still 
present,  steps  (b)  and  (c)  are  repeated  until  no 
appreciable  difference  in  magnitude  and  distribution 
of  stresses  is  obtained  upon  further  iteration.  It 
is  important  to  recognize  that  the  initial  elastic 
stiffness  matrix  is  used  throughout  the  computation. 

This  type  of  analysis  has  been  used  by  Zienkiewicz,  et  al .  (1968) 
in  the  solution  of  plane  strain  problems  in  rock  mechanics.  Heuze, 
Goodman,  and  Bornstein  (1971)  have  used  the  same  technique  to  study 
the  deformability  of  jointed  rock  in  borehole  jacking  tests.  The 
convergence  of  the  solution  has  been  found  to  be  very  slow.  The 
number  of  iterations  required  for  convergence  ranged  from  10  to 
15  cycles.  Zienkiewicz,  et  al .  (1968)  reported  that  faster  con¬ 
vergence  could  be  obtained  if  the  elastic  constants  are  modified 
by  reducing  the  modulus  in  the  direction  of  the  principal  tensile 
stress.  With  this  modification,  the  stiffness  matrix  would  have 
to  be  recomputed  at  every  stage. 

2 .  "Joint  Perturbation"  Analysis 

A  one-dimensional  joint  element  was  developed  by  Goodman,  Taylor, 
and  Brekke  (1968)  and  applied  to  several  rock  mechanics  problems. 
This  technique  is  capable  of  modelling  the  behavior  of  joints, 
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bedding  planes  and  other  geologic  discontinuities, 
version  defined  normal  and  shear  stiffnesses,  and  Kg,  for 
joint  elements  and  incorporated  the  stiffness  of  these  elements 
into  the  stiffness  of  the  overall  structure.  Fig.  1  shows  a 
one-dimensional  joint  element  with  its  local  coordinate  system 
and  the  sign  convention  on  stresses  and  displacements.  The 
joint  element  was  assumed  to  have  no  thickness.  The  original 
version  was  only  valid  for  a  linear  elastic  analysis.  In  his 
study  on  the  deformability  of  joints,  Goodman  (1969)  has  reported 
that  there  are  four  typical  shear  stress -deformation  relationships 
for  various  weak  surfaces  as  shown  in  Fig.  2.  It  is  apparent  that 
peak  shear  strength  in  most  cases  is  greater  than  residual  shear 
strength.  After  the  peak  shear  strength  is  reached,  shear  stress 
drops  to  residual  values  and  appreciable  movement  can  take  place 
without  increase  in  shear  stress.  It  is  obvious  that  this  type 
of  stress -deformation  behavior  cannot  be  approximated  by  one 
single  value  of  shear  stiffness  as  used  in  the  original  analysis, 
Goodman  et  al .  (1968).  Heuze,  Goodman  and  Bornstein  (1971)  modi¬ 
fied  the  original  version  so  that  normal  and  tangential  properties 
can  be  varied  as  a  function  of  displacement.  They  used  an  itera¬ 
tive  approach  to  solve  this  type  of  nonlinear  problem.  In  their 
procedure,  the  same  boundary  value  problem  is  analyzed  repeatedly. 
For  each  analysis,  a  new  set  of  normal  and  tangential  stresses 
and  displacements  across  a  joint  element  is  calculated  and  used 
together  with  joint  properties  and  joint  strength  to  modify  the 
normal  and  shear  stiffness  for  the  next  iteration.  The  itera¬ 
tive  technique  used  by  Heuze,  Goodman  and  Bornstein  (1971)  is 
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FIG  1  -  LINKAGE  OR  "JOINT"  ELEMENT  WITH 
ITS  LOCAL  COORDINATE  SYSTEM 
(After  Goodman,  Taylor  and  Brekke ,  1968) 
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FIG.  2  -  TYPICAL  SHEAR  STRESS  -  DEFORMATION  RELATIONSHIPS  FOR  VARIOUS 
WEAKNESS  SURFACES  (AFTER  GOODMAN,  1969) 


illustrated  in  Fig.  3.  In  this  approach,  a  new  stiffness  matrix 
has  to  be  formulated  and  the  same  boundary  value  problem  is  solved 
for  every  iteration.  This  is  a  significant  computational  differ 
ence  from  the  stress  transfer  technique  used  in  the  no  tension 
analysis  of  Zienkiewicz  et  al.  (1968)  where  the  stiffness  matrix 
is  not  changed  with  repeated  iterations. 

Zienkiewicz,  et  al.  (1970)  have  shown  that  an  iterative  process 
similar  to  "stress  transfer"  analysis  may  be  employed  in  the 
analysis  of  jointed  rock  systems. 

3.  Elasto-plastic  Analysis 

An  elasto-plastic  analysis  has  been  suggested  by  various  re¬ 
searchers  to  account  for  the  possible  yielding  of  rock  due  to 
the  stress  concentrations  induced  in  the  rock  mass  around  an 
underground  opening.  Reyes  and  Deere  (1966)  developed  a  method 
based  on  the  incremental  theory  of  plasticity  to  study  elasto- 
plastic  behavior  of  underground  openings.  From  the  computational 
point  of  view,  the  process  used  by  Reyes  and  Deere  (1966)  has  the 
disadvantage  in  that,  at  each  iteration  the  stiffness  of  the 
structure  is  changed,  requiring  a  reformulation  and  computation 
of  the  stiffness  matrix  which  will  involve  extra  computer  time. 

An  alternative  approach  has  been  developed  by  Zienkiewicz,  et  al. 
(1969).  This  approach  uses  a  technique  referred  to  as  the  "ini¬ 
tial-stress"  technique  which  is  consistent  with  the  "no  tension" 
analysis.  The  solution  obtained  by  the  "initial  stress"  process 
involves  a  series  of  load  increments.  For  each  load  increment, 
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FIG.  3  -  ITERATIVE  TECHNIQUE  FOR  JOINT  PERTURBATION  ANALYSIS 
(AFTER  HEUZE,  GOODMAN  AND  BORNSTEIN,  1971) 
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the  solution  satisfying  equilibrium  and  yield  criteria  is 
achieved  by  a  series  of  approximations  or  iterations.  Fig.  4 
illustrates  how  the  final  solution  may  be  achieved  by  a  series 
of  iterations.  Baker,  Sandhu,  and  Shieh  (1969)  have  also 
developed  a  computational  technique  fcor  the  elasto -plastic 
analysis.  The  basic  concept  is  quite  similar  to  that  used  in 
the  initial  stress  approach  presented  by  Zienkiewicz,  et  al. 
(1969).  To  save  computer  execution  time,  both  approaches  use 
the  initial  (elastic)  stiffness  of  the  structure  at  each  step 
of  computation.  Corrective  body  forces  are  calculated  and 
applied  to  the  structure  during  iterations  to  insure  that  the 
element  is  just  on  the  yield  surface.  The  only  difference 
between  the  approaches  of  Zienkiewicz,  et  al .  and  Baker,  et  al. 
is  that  the  latter  employed  a  different  method  to  insure  the 
convergence  of  the  final  solution  for  each  load  increment. 

4 .  Time -dependent  Analysis 

Time -dependent  analysis  for  problems  in  rock  mechanics  can  be 
considered  in  two  categories.  The  first  includes  those  cases 
where  the  boundary  value  problem  changes  with  time.  Such  prob¬ 
lems  include  consideration  of  the  gradual  (time -dependent) 
creation  of  the  excavation  or  the  installation  of  reinforcement 
(e.g.,  rock  bolts)  at  various  times  during  the  construction  of 
the  opening  or  later  during  its  performance.  The  second  cate¬ 
gory  includes  those  problems  where  the  properties  of  the  rock 
surrounding  the  excavation  are  time-dependent  (creep) .  This 
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FIG.  4  -  GRAPHICAL  INTERPRETATION  OF  THE  II 
PROCESS  FOR  ELASTO- PLASTIC  ANALYS 
ZIENKIEWICZ ,  ET  AL  1969) 
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review  was  confined  to  problems  in  the  second  category,  and  in 
the  subsequent  discussions,  time -dependent  analyses  refer  only 
to  time-dependent  material  properties. 

Not  much  emphasis  has  been  placed  on  developing  time -dependent 
analyses  for  rock  mechanics  problems  because  the  great  majority 
of  rocks  do  not  exhibit  significant  time -dependent  behavior. 

In  structural  analysis,  two  approaches  have  been  taken  in  devel¬ 
oping  solutions  for  time -dependent  problems.  These  depend  on 
the  complexity  of  material  response.  For  linear  problems,  the 
theory  of  linear  viscoelasticity  has  been  utilized.  For  mate¬ 
rials  exhibiting  nonlinear  behavior,  empirical  stress -strain 
relations  together  with  iterative  solution  techniques  have  been 
used.  Excessive  computer  costs  and  difficulties  associated  with 
determining  reliable  material  properties  have  limited  the  appli¬ 
cation  of  such  analytical  techniques  to  problems  in  rock  mechanics. 

Time-dependent  analyses  for  nonlinear  material  properties  have 
been  developed  by  Deere  and  Boresi  (1963),  Nair  (1967),  Aiyer 
(1969),  and  Nair  and  Boresi  (1970).  All  these  analyses  were 
developed  for  specialized  problems  and  with  the  exception  of 
Aiyer's  work  were  confined  to  spherically  symmetric  or  axisym- 
metric  problems.  However,  at  the  time  of  this  review,  no  pub¬ 
lished  information  was  found  on  the  availability  of  a  plane 
finite  element  program  which  had  been  utilized  to  analyze  prob¬ 
lems  in  rock  mechanics.  The  approach  utilized  by  Nair  and 
Boresi  (1970)  uses  the  finite  element  method  and  uses  a  computa¬ 
tional  technique  which  can  be  applied  to  the  analysis  of  plane 
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problems.  The  basic  concept  used  in  this  analysis  is  similar 
to  that  proposed  by  Greenbaum  (1966)  .  The  method  of  incorpo¬ 
rating  the  time -dependent  behavior  of  the  material  is  briefly 
described  as  follows:  As  a  first  step  in  the  computation,  an 
elastic  stress  distribution  is  computed;  the  effective  stress 
is  then  computed.  Based  on  the  value  of  the  effective  stress, 
the  effective  strain  rate  is  computed  using  the  creep  stress- 
strain  time  law.  Taking  a  small  time  increment  At,  the  incre¬ 
mental  strains  can  now  be  computed.  From  these  incremental 
strains,  the  incremental  stresses  are  computed  using  the  elastic 
stress -strain  relations.  These  incremental  stresses  are  then 
input  into  the  program  as  initial  stresses,  converted  into  nodal 
point  forces,  and  the  problem  is  solved  as  an  elastic  problem 
and  the  new  stresses  and  displacements  are  determined.  These 
stresses  are  then  used  to  develop  the  new  incremental  strains 
and  the  process  repeated.  In  this  manner,  the  variation  of 
stress  and  displacement  with  time  is  determined. 

Summary  of  Achievement  for  Phase  la 

A  review  of  the  available  literature  indicated  that  methods  are 
available  for  modelling  geologic  discontinuities  and  failure 
characteristics  that  occur  in  rock  masses.  However,  the  com¬ 
putational  techniques  associated  with  the  various  methods  are 
not  identical.  Consequently,  the  first  step  in  the  development 
of  any  general  purpose  computer  program,  which  would  include  the 
capabilities  of  the  individual  programs,  would  be  to  make  the 
various  computational  techniques  consistent. 
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It  was  also  found  that  time -dependent  analyses  were  not  developed 
to  the  level  of  sophistication  of  the  time - independent  analyses. 
Furthermore,  the  cost  associated  with  a  time -dependent  analysis 
and  the  fact  that  relatively  few  rocks  exhibit  significant  time- 
dependent  behavior  leads  to  the  conclusion,  at  this  time,  that 
the  development  of  a  general  purpose  program  for  the  time- 
independent  plane  problems  be  treated  separately  from  the  devel¬ 
opment  of  a  program  for  time -dependent  analyses.  In  this  study, 
only  the  former  has  been  considered. 

DEVELOPMENT  OF  A  GENERAL  COMPUTER  PROGRAM  (Phase  lb) 

In  order  to  develop  a  single  general  computer  program,  it  was 
necessary  to  develop  a  consistent  computational  technique  to 
model  the  different  aspects  of  rock  behavior.  On  the  basis  of 
the  review  of  the  available  techniques,  it  was  concluded  that 
the  "initial  stress"  (stress  transfer)  technique  presented  by 
Zienkiewicz  and  his  co-workers  would  provide  a  consistent  approach 
in  the  development  of  a  consolidated  computer  program,  for  the 
stress  analysis  of  plane  problems,  which  is  capable  of  modelling 
joints  and  other  geologic  discontinuities,  and  elasto -plastic  or 
time-dependent  rock  properties.  The  major  reason  for  selecting 
this  approach  was  the  computational  advantage  that  results  from 
using  the  initial  elastic  stiffness  at  every  stage  of  the  solu¬ 
tion  process.  The  features  incorporated  in  the  combined  computer 
program  for  time - independent  plane  problems  are:  (1)  No  Tension 
analysis,  (2)  Joint  Perturbation  analysis,  and  (3)  Elasto-plastic 
analysis.  In  order  to  develop  a  single  computer  program  using 
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the  selected  computational  technique,  it  was  necessary  to 
formulate,  write  and  debug  programs  for  joint  perturbation 
and  elasto-plastic  analysis.  The  essential  concepts  used  to 
include  the  above  listed  rock  characteristics  are  discussed  in 
the  subsequent  paragraphs. 

I .  No  Tension  Analysis 

The  basic  concept  used  in  the  combined  computer  program  for 
the  no  tension  analysis  is  similar  to  that  developed  by 
Zienkiewicz,  et  al.  (1968).  The  major  steps  in  the  analysis 
used  can  be  summarized  as  follows: 

1.  Assign  initial  stresses  to  the  rock  mass,  and  cal¬ 
culate  the  boundary  loads  required  on  the  cavity 
face  to  simulate  the  creation  of  the  opening  and 
other  structural  loads  applied  to  the  system. 

2.  Analyze  the  problem  as  a  linear  elastic  problem. 

Add  the  induced  changes  in  stress  to  the  initial 
stresses  and  compute  the  principal  stresses. 

3.  Determine  those  elements  in  which  tension  exists. 

If  the  material  is  assumed  incapable  of  sustaining 
tension,  or  if  the  tensile  stress  exceeds  the  tensile 
strength,  it  is  necessary  to  eliminate  the  excess 
tensile  stresses.  This  is  done  by  applying  nodal 
point  forces  calculated  to  eliminate  the  excess 
tensile  stresses. 
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4.  The  elastic  analysis  is  repeated  for  the  calculated 
equivalent  nodal  point  forces;  stresses  are  determined 
in  the  elements.  The  check  for  tensile  stresses  is 
repeated . 

5.  If,  at  the  end  of  step  (4),  principal  tensile  stresses 
are  still  present,  steps  (3)  and  (4)  are  repeated  until 
there  is  no  appreciable  difference  in  magnitude  and 
distribution  of  stresses  with  further  iterations. 

It  has  been  reported  by  Nair  and  Chang  (1971),  and  Sandhu ,  et  al . 
(1971)  that  using  the  technique  proposed  by  Zienkiewicz,  et  al . 
(1968)  ,  a  large  number  of  iterations  is  required  to  reduce  ten¬ 
sile  stresses  to  negligible  values.  Furthermore,  the  solution 
may  not  converge  in  some  cases  where  one  principal  stress  is 
tensile  and  the  other  compressive.  The  following  modifications 
were  made  in  the  analysis  to  increase  the  rate  of  convergence. 

In  Step  (3) ,  when  the  linear  elastic  solution  indicates  that 
the  rock  is  subjected  to  tensile  stresses  greater  than  the 
tensile  strength,  the  rock  is  assumed  fractured  and  incapable 
of  transferring  stresses  between  two  orthogonal  directions.  In 
transferring  tensile  stresses,  the  Poisson's  ratio  for  that 
element  is  reduced  to  zero.  This  technique  has  been  employed 
by  Nair  and  Chang  (1971)  in  conjunction  with  an  iterative  pro¬ 
cess  where  the  stiffness  matrix  was  modified  at  each  state.  In 
the  present  study,  since  the  same  stiffness  matrix  is  used  in 
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all  iterations,  to  consider  the  Poisson's  effect  a  correction 
is  made  on  the  stress  before  the  stress  transfer  process  is 
performed.  As  shown  in  Fig.  5,  where  only  the  tensile  stress 
in  one  direction  (a  ')  is  to  be  transferred,  a  stress  state  in 
accordance  with  the  following  equations  has  to  be  applied  in 
order  to  eliminate  Poisson's  effect. 


a  =  a  '  +  v2a  '  +  .  .  . 

A  A  A 

Oy  -  vax'  +  v3ax'  +  .  .  . 


CD 


In  the  case  where  both  principal  tensile  stresses  ox'  and  o' 
are  to  be  transferred,  the  state  of  stress 


ox  =  (1  +  v2  +  . ..)  ox  2  +  v  (1  +  v2  +  ...)  o'y2 
Oy  =  (1  +  v2  +  ...)  Oy 2  +  v  (1  +  v2  +  ...)  o'x2 


(2) 


is  applied  to  the  element.  The  rate  of  convergence  is  increased 
when  this  computational  procedure  is  utilized. 


II .  Joint  Perturbation  Analysis 

The  stiffness  matrix  of  one -dimensional  joints  is  formulated 
according  to  the  procedure  developed  by  Goodman ,  Taylor,  and 
Brekke  (1968).  The  one -dimensional  joint  with  its  local  coor¬ 
dinate  system  and  the  sign  convention  of  its  normal  and  shear 
stresses  have  been  presented  in  Fig.  1.  The  approach  used  in 
this  study  for  the  nonlinear  analysis  of  joint  elements  is 
similar  to  the  stress  transfer  technique  proposed  by  Zienkiewicz 
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FIG.  5  -  STRESS  TRANSFER  TECHNIQUE  TO  ELIMINATE  POISSON'S 
EFFECT  FOR  ELEMENTS  IN  TENSION 
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and  his  co-workers.  The  shear  strength  of  a  joint  element  depends 
on  the  cohesion,  c,  the  friction  angle,  <j>,  the  joint  roughness  and 
the  normal  stress  acting  across  the  joint.  Patton  (1966)  has  shown 
that  the  influence  of  joint  roughness  on  the  shear  strength  along 
a  rock  surface  -can  be  taken  into  account  by  increasing  the  friction 
angle,  <j>,  of  the  joint  surface  by  an  amount  <j>r,  which  is  the  average 
angle  between  the  undulations  on  the  joint  surface  and  the  direction 
of  sliding  along  the  joint.  Therefore,  the  effective  friction  angle 
<p  of  a  rough  joint  surface  is  given  by 

v 

<l>e  =  <f>  +  <f>r  (3) 

and  the  shear  strength  of  a  joint  may  be  expressed  by 

t£  =  o  +  aN  tan  (4) 

where 

t£  =  shear  strength  of  the  joint 
c  =  cohesion  along  the  joint 

*  normal  stress  across  the  joint 
<f>  =  effective  friction  angle  of  the  joint  surface. 

If  the  normal  stress  across  the  joint  is  tensile,  it  is  assumed 
that  the  joint  is  incapable  of  resisting  any  shear  stress,  i.e., 
it  has  no  strength. 

The  procedure  used  to  account  for  the  nonlinear  behavior  of  the 
joint  elements  is  as  follows: 
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1.  Initial  stresses  are  assigned  to  all  elements,  and 
the  boundary  loads  required  on  the  cavity  face  to 
simulate  the  creation  of  the  opening  and  other  loads 
applied  to  the  system  are  calculated. 

2.  The  problem  is  analyzed  as  a  linear  elastic  problem 
From  the  computed  nodal  point  displacements,  the 
normal  and  tangential  displacements  across  the  joint 
are  calculated.  These  are  used  to  compute  the  changes 
in  normal  and  tangential  stresses  using  the  normal  and 
shear  stiffness  of  the  joint.  The  induced  changes  in 
normal  and  tangential  stresses  are  combined  with  the 
initial  stresses. 

3a.  Since  the  joint  is  assumed  incapable  of  sustaining  any 
shear  stress  under  a  tensile  normal  stress,  a  check  is 
made  to  see  if  a  tensile  normal  stress  exists  across  the 
joint.  If  a  tensile  normal  stress  is  present,  then  both 
normal  and  tangential  stresses  are  eliminated  and  the 
equivalent  nodal  point  forces  around  the  joint  are 
calculated . 

3b.  If  the  normal  stress  is  compressive,  the  shear  strength 
of  the  joint  is  calculated  by  Equation  (4) .  If  the  tan¬ 
gential  stress  is  less  than  the  calculated  shear  strength, 
the  joint  remains  intact.  The  next  step  is  to  check  clo¬ 
sure  as  in  step  3(c).  If  the  tangential  stress  is  greater 
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than  the  calculated  shear  strength,  the  joint  is  in 
yield.  The  excess  tangential  stress  which  is  the 
difference  between  the  tangential  stress  and  shear 
strength  is  eliminated  and  replaced  by  equivalent 
nodal  point  forces. 

3c.  If  the  joint  is  in  compression,  and  it  undergoes  excess 
closure  beyond  the  allowable  closure  specified,  the 
nodal  point  displacements  around  the  joint  are  corrected 
to  limit  the  displacement  to  the  maximum  allowable  clo¬ 
sure.  Equivalent  normal  forces  for  the  corrections  in 
the  nodal  point  displacements  are  computed. 

4.  The  elastic  analysis  is  repeated  for  the  equivalent 
nodal  point  forces  calculated  in  Step  (3).  The  initial 
stiffness  of  the  system  is  used  throughout  the  computa¬ 
tion.  On  the  basis  of  the  computed  stresses  and  dis¬ 
placements,  the  checks  described  in  Step  (3)  are  repeated. 

5.  If,  on  repeated  Step  (3),  corrective  nodal  point  forces 
around  the  joints  are  still  present,  the  analysis  pro¬ 
ceeds  to  Step  (4).  This  iterative  process  is  continued 
until  the  change  in  magnitude  and  distribution  of  stresses 
obtained  upon  iteration  is  negligible. 

Fig.  6  illustrates  schematically  the  stress  transfer  process  for 
the  nonlinear  analysis  of  the  joint  systems  described  above. 
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A  flow  diagram  illustrating  the  stress  transfer  technique  for 
joint  perturbation  analysis  is  shown  in  Fig.  7.  It  should  be 
noted  that  in  this  technique,  the  joint  stiffnesses  are  held 
constant  while  transferring  excess  stresses  across  joint  elements. 
Heuze,  Goodman  and  Bornstein  (1971)  use  a  different  technique  in 
conducting  the  joint  perturbation  analysis  in  that  they  employ 
an  iterative  process  in  which  the  joint  stiffnesses  are  adjusted 
in  each  iteration. 

HI.  Elasto -Plastic  Analysis 

In  developing  an  elasto -plastic  analysis,  it  is  necessary  to 
define  a  yield  function  and  the  stress -strain  relations  before 
and  after  yield.  Prior  to  yield,  it  is  assumed  that  linear 
elastic  stress -stra’ n  relations  are  applicable. 

Yield  Function 

in  the  present  study,  the  yield  function  utilized  is  a  generaliza¬ 
tion  of  the  Mohr -Coulomb  hypothesis  suggested  by  Drucker  and  Prager 
(1952).  The  yield  function  is  represented  by  the  following  equation 

f  =  aJj  +  /J7  =  k  (5) 

where : 

l 

ot  and  k  =  material  constants 

Ji  =  first  stress  invariant 
J  2  =  second  invariant  of  stress  deviation 

The  stress  invariants  may  be  expressed  in  terms  of  the  stress 
components  as  follows: 

J.  =  °x  +  °v  *  °z  (6) 
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Initiat  stresses  ca/cu  lated  and  boundary  loads  applied 


So  ire  /-ha.  boundary  ra/ue.  prob  /««  a-s 
a.  //near  e /as  tic  case. 


Compute  normal  and  tangential  displace¬ 
ment  and  stresses  for  Joint  elements 


A _ _ /  rt>»crl  Compressi V 


J4f.r  I  Ca.lcula.ie  shear  strength,  Xf 
If  tangential  stress  X*  ~^ft 


Calculate  excess  tangential 
stress  to  be  redistributed 

z-tf 


Compute  e<ju  /  t/brating 
Poda. f  point  -forces 


Jf  the  e /ement  undergoing 
closure  /'. e .  norma/  displacement, 
!/>  max.  a//ocu<xb/e,  closure ,  I/max.? 


Redistribute  excess  stresses 


FIG,  7  -"STRESS  TRANSFER"  TECHNIQUE  FOR  JOINT  PERTURBATION  ANALYSIS 
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J2  =  t(ox  •  V  +  (oy  '  0z)2  +  (oz  '  °x)2]  +  T 


xy 


2  2 
+  t  +  T 
y  z  zx 


(7) 


The  yield  surface  expressed  by  Equation  (5)  for  a  >  0  is  a 
right  circular  cone  with  its  axis  equally  inclined  to  the 
coordinate  axes.  For  a  =  0  equation  (5)  reduces  to  the  Von 
Mises  yield  function;  the  yield  surface  is  a  right  circular 
cylinder.  These  two  yield  surfaces  are  illustrated  in  Figure  8. 


In  the  case  of  plane  strain,  Drucker  and  Prager  (1952)  have 
shown  that 


a 


_ tan  4> _ 

(9  +  12  tan2  4>) ^ 


and 


(9  +  12  tan2  <j>) ^ 

where : 

c  3  the  cohesion  of  the  material 

<}>  3  angle  of  internal  friction  of  the  material 


(8) 


(9) 


Stress-Strain  Relations 

During  an  infinitesimal  increment  of  stress,  changes  in  strain 
are  assumed  to  be  composed  of  elastic  and  plastic  parts  if  the 
element  is  in  yield,  i.e. 


{Ae}  3  {Ae}  +  {Ae} 

e  p 


(10) 
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FIG.  8  -  GENERALIZED  MOHR- COULOMB  YIELD  SURFACE 
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The  elastic  strain  increments  are  related  to  stress  increments 
by  the  generalized  Hooke’s  law  as 


{Ao}  *  [D]  {Ae} 


(ID 


where  the  strain-stress  matrix  is 


£D?  =  TTVv)  (1  -  2\jy 


(1  -  V)  V  0 
V  (1  -  V)  o 
0 


o  (1  -  2v) 

2 - Jj 


(12) 


and  E  is  the  elastic  modulus  and  v>  the  Poisson's  ratio  for  the 
linear  isotropic  elastic  material. 

Utilizing  the  Drucker  Prager  criteria  (equation  5)  Reyes  (1966) 
developed  elasto-plastic  stress-strain  relations.  These  relations 
can  be  expressed  as  follows: 


(AO)  =  [D]e.p.Ue)p 


(13) 


where 


CD] 


e.p 


'll 

D12 

D13 

'21 

D22 

D23 

31 

D32 

D33 
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D  =  2 G ( 1  -  h  -  2h  a  -  h  a  2) 

11  2  1  X  3  X 


D  =  2G (1  -  h  -  2h  a  -  h  a,2) 
22  2  i  y  3  y 


3  3 


=  2G(I  ■  h 3  Txy2) 


(14) 


D  =  D  =  -2G  [h  +  h  (a  +  av)  +  h  a  a  ] 
12  2i  2  ixy  3xy 


D  =  D  =  -2G(h  t  +  h  a  t  ) 
13  3i  i  xy  3  x  xy 


D  =  D  =  -2G  (h  t  +  h  a  t  ) 

2  3  3  2  1  xy  3  y  xy 


where : 


B 


G  (shear  modulus)  =  2~(l'~+"  v) 


h  = 

i 


3Ka 

xr 


,2  K> 


(15) 


h  = 

2 


J 

3Ka  Ji 

6  J1 

G  3  J  ^ 

_  2  _ 

„  2  J 

3  v  K  k 


(1  +  9a2  §  ) 


E  J  55  (1  +  9a2  £) 


h 

3 


_ 1 

2  J  (1  +  9a2  ^  ) 


K  (bulk  modulus) 


E 

=  3(1  -  2vJ 


Since  the  elasto-plastic  stress-strain  relation  is  a  function  of 
the  elastic  constants,  the  yield  parameters  (a,  k)  and  the  stress 
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state,  it  is  necessary  to  keep  a  record  of  the  axial  stress  a  . 

For  the  plane  strain  case  the  change  in  stress, Ac  ,  is  related 

Ld 

to  the  changes  in  stress  in  x  -  y  plane  in  the  elastic  range  by 

A cz  =  v(Aa^  (16) 

where  Ac^  and  Ac^  are  changes  in  stress  in  two  principal  stress 
directions  in  x  -  y  plane.  In  the  plastic  range,  Drucker  and 
Prager  (1952)  have  shown  that  the  change  in  axial  stress  is 
given  by 

A az  =  j  (A +  Ac^)  -  j  (Aa^  -  Aa ^ )  Sin  <j>  (17) 

Computational  Procedure 

To  conduct  an  incremental  elasto-plastic  analysis  the  initial 
stress  approach  developed  by  Zienkiewicz,  et  al  (1969)  is  utilized. 
This  approach  is  consistent  with  that  described  previously  for  the 
no  tension  and  joint  perturbation  analyses.  The  load  is  applied  in 
a  series  of  increments.  The  initial  stress  process  approaches  the 
solution  of  a  nonlinear  problem  by  a  series  of  approximations. 

The  procedure  for  conducting  the  elasto-plastic  analysis  during 
a  typical  load  increment  can  be  summarized  as  follows: 

1.  For  the  applied  load  increment,  the  elastic  increments  of 

stress  { Aa ' } ^  and  corresponding  strain  {Ae'}^  are  determined. 
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2.  Add  {Aa1}-^  to  the  stresses  existing  at  the  end  of  the  last 

increment  {aQ}  to  obtain  {a1}.  The  corresponding  change  in 

axial  stress  Act  '  is  calculated  according  to  one  of  the  fol- 
z 

lowing  two  conditions. 

(a)  If  f(aQ)  -  k  ^  0  i.e.,  the  element  was  in  yield  at  the 
start  of  increment 

A  o’  =  j  .( Aa  +  Aa2)  -  j  (Aa^  -  Ac?2)  sin  (j>  (18) 

(b)  If  f  (aQ)  -  k  <  0  i.e.,  the  element  was  in  the  elastic 
range  at  the  start  of  increment 

A  o'  =  v  (  Aa^  +  Aa2)  (19) 

and  the  current  axial  stress  is  obtained  by 


a 


! 

z 


Caz)0  + 


Aa  ' 
z 


(20) 


Calculate  [f(a')  -  k]  using  equations  (5)  through  (9).  If 
f(aQ)  -  k  <  0  and  f(a')  -  k  <  0,  only  changes  in  elastic  strain 
occur  (i.e.,  there  is  no  yield).  No  further  computations 
are  required.  If  f(cQ)  -  k  >_  0  and  f(c')  -  k  <  0,  this 
implies  that  no  further  yielding  occurred  during  the  applied 
load  increment.  The  corresponding  change  in  axial  stress  is 
corrected  using  equation  (19)  and  a  '  recalculated  by 
equation  (20) 
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3.  I  £  f(a')  -  k  0  and  f(aQ)  -  lc  =  0  i.e.,  the  element  was 
in  yield  at  the  start  of  increment,  calculate  {Aa}^,  using 
the  following  relation: 

l4o}l  ■  tDle.p.(4£,'l  (21) 

where  [D]  is  the  elasto -plastic  stress -strain  relation 

c  •  p  • 

expressed  by  equation  (14) . 

The  excess  stresses  which  have  to  be  redistributed  are 
calculated  as  follows: 

{Aa"}^  =  {Aa'J^  -  {Aa}.^  (22) 

whereas  (Aa  ")]_  is  computed  in  accordance  with  the  following 

A 

equation : 

CAaz")i  =  j  (Aa-^"  +  A02")  -  j  (Aa^"  -  Ac^")  sin  4  (^) 

The  current  stresses,  which  are  stored,  are  computed  as 
follows : 

(a }  =  (o'  }  -  {AaM)1  (24) 

and  current  strains 

{e}  =  {eq}  +  (Ae ' (25) 
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4.  If  f  (0 ' )  -  k  >  0  but  £ (aQ)  -  k  <  0  i.e.,  the  element  goes 

from  the  elastic  to  the  plastic  range,  and  it  is  necessary  to 
find  the  intermediate  stress  value  at  which  yielding  commences. 
This  is  done  by  the  following  interpolation  procedure: 


Af  -fCUo'Jj)  (26) 

{ Ao  ' } ,  o  „  =  A  {Ao'k  C27) 

l  0  •  P  •  -k 


where 


A  =  -  k 

A  aT 


{0'}  =  (o’ } 

e .  p . 


{Ao')1 


e .  p . 


(28) 


Use  {0'}  e.p.  t0  compute  [ D ]  g  ^  and  equation  (21)  to  calcu¬ 
late  { Ao}^ .  Then  proceed  as  in  step  (3). 


5.  Considering  {Ao"}-^  as  initial  stresses,  the  corresponding 
equilibrating  nodal  point  forces  {P}e^  are  computed. 


6.  The  system  is  solved  for  the  loads  {P}^  and  { A0 ' } 2  and 
{Ac1 >2  are  determined. 

7.  Steps  2  to  6  are  repeated  until  the  nodal  forces  {P}  reach 
sufficiently  small  values. 

A  flow  diagram  showing  the  stress  transfer  technique  for  elasto- 
plastic  analysis  is  illustrated  in  Figure  9. 
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FIG.  9  -  "STRESS  TRANSFER"  TECHNIQUE  FOR  ELASTO-PLASTIC  ANALYSIS 
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Experience  from  a  limited  number  of  solutions  performed  indicated 
that  the  rate  of  solution  convergence  would  depend  on  the  stress 
state  and  the  yield  constants  a  and  k  given.  In  general,  it  was 
found  that  if  the  Von  Mises  yield  criterion  is  used  i.e.,  a  =  0 
or  <j>  =  0,  the  solution  convergence  is  rapid,  three  to  six  itera¬ 
tions  being  necessary  in  any  increment.  However,  when  <|>  is  non¬ 
zero  the  solution  converges  considerably  slower.  This  convergence 
problem  was  also  noted  by  Baker,  Sandhu  and  Shieh  (1969).  In  some 
cases,  it  may  be  found  necessary  to  use  various  numerical  schemes 
to  accelerate  convergence. 

Summary  of  Achievement  for  Phase  lb 

A  computational  procedure  using  the  initial  stress  (stress  transfer) 
approach  developed  by  Zienkiewicz  and  his  co-workers  was  formulated 
for  modelling  certain  classes  of  joints,  no  tension  characteristics 
in  the  rock  mass  and  elasto-plast ic  behavior.  This  procedure  was 
used  to  develop  a  computer  program  which  could  determine  the  mechan¬ 
ical  state  in  a  rock  mass  with  the  above-mentioned  characteristics 
in  the  vicinity  of  an  excavation.  Because  available  programs  did 
not  use  a  single  computational  technique,  it  was  necessary  to  write 
new  programs . 

ILLUSTRATIVE  PROBLEMS  (Phase  lc) 

Definition  of  Problems 

The  following  examples  were  analyzed  using  the  combined  finite 
element  computer  program  developed  in  this  study  for  the  purpose 
of  verifying  and  demonstrating  the  use  of  the  program. 
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I.  Elasto-plastic  analysis  of  a  thick-walled  circular  tube 
with  the  Von  Mises  yield  criterion.  A  closed  form  solu¬ 
tion  is  available  for  this  case  for  verification. 

II.  Elasto-plastic  analysis  of  a  circular  opening  with  the 
generalized  Mohr-Coulomb  yield  criterion.  The  results 
are  compared  with  those  obtained  by  Reyes  (1966)  and 
Baker,  et  al  (1969) , 

III.  Combined  no  tension,  joint  perturbation  and  elasto-plastic 
analysis  of  a  rectangular  underground  opening  to  demon¬ 
strate  the  usage  of  the  combined  computer  program. 

Results 

I.  Elasto-Plastic  Analysis  of  a  Thick-walled  Circular  Tube 
with  the  Von  Mises  Yield  Criterion  Subject  to  Internal 
Pressure 

Prager  and  Hodge  (1951)  have  presented  a  closed  form  solution  of 
the  above  problem.  The  material  is  assumed  to  obey  the  Von  Mises 
yield  criterion,  a  special  case  of  the  Mohr-Coulomb  criterion  ex¬ 
pressed  by  Eq.  (5)  in  which  <J>  =  0  or  a  =  0  and  k  =  c.  The  dimen¬ 

sions  of  the  tube,  the  material  properties  and  the  finite  element 
idealization  of  the  problem  are  shown  in  Figure  10. 

The  circular  tube  was  subjected  to  internal  pressure  up  to 
p  *  1,39  k  in  five  unequal  increments.  The  results  of  the 
analysis  indicated  that  the  interior  surface  reached  the  yield 
surface  after  the  first  increment  of  loading  (p  >  0.75k). 
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FIG.  10  -  FINITE  ELEMENT  MESH  FOR  AN  ELASTO- PLASTIC  ANALYSIS  OF  A 
THICK-WALLED  CIRCULAR  TUBE  (b  =  2a) 


-39- 


Subsequent  increase  in  pressure  caused  the  material  near  the 
interior  surface  to  undergo  elasto-plastic  behavior.  For  each 
subsequent  pressure  increment,  six  cycles  of  stress  redistribution 
were  required  for  solution  convergence,  except  for  the  second  in¬ 
crement  for  which  only  three  cycles  were  required.  The  results 
of  the  analysis,  together  with  the  closed  form  solution  are  shown 
in  the  form  of  stresses  and  displacements  in  Figs.  11  through  14. 
Fig.  11  presents  the  distribution  of  circumferential  stress  for 
the  various  pressure  increments.  The  apexes  of  these  computed  curves 
indicate  the  boundary  of  the  plastic  zone.  It  should  be  recognized 
that  the  values  of  stress  are  plotted  at  the  centroid  of  the  element 
From  Fig.  11,  the  ratio  of  the  radius  of  the  plastic  zone  boundary 
to  the  internal  radius  (p/a)  is  computed  for  each  pressure  increment 

Figs.  12  and  13  present  the  distribution  of  radial  and  axial  stress 
for  all  pressure  increments.  Also  shown  on  these  figures  is  the 
p/a  corresponding  to  the  distribution.  The  results  can  be  compared 
with  the  closed  form  solution  for  the  corresponding  p/a.  The 
results  are  in  excellent  agreement.  Fig.  14  shows  the  radial  dis¬ 
placements  on  the  interior  surface  as  functions  of  the  radius  p  of 
the  elastic-plastic  boundary.  Comparison  between  the  results 
obtained  from  the  finite  element  analysis  and  those  from  the 
closed  form  solution  indicates  good  agreement. 

II.  Elasto-plastic  Analysis  of  a  Circular  Opening  with  the 
Generalized  Mohr-Coulomb  Yield  Criterion 
For  purposes  of  demonstrating  an  elasto-plastic  analysis  for  a 
material  with  the  generalized  Mohr-Coulomb  yield  criterion,  a 
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FIG.  11  -  DISTRIBUTION  OF  CIRCUMFERENTIAL  STRESS 
FOR  A  THICK -WALLED  CIRCULAR  TUBE 


•  -  INCREMENT  1  p 

O  -  INCREMENT  2,  ITERATION  3,  p 

A  -  INCREMENT  3,  ITERATION  6,  p 

□  -  INCREMENT  4,  ITERATION  6,  p 

0  -  INCREMENT  5,  ITERATION  6, 

p  =  1.39k 

p  =  RADIUS  TO  ELASTO- 
PLASTIC  BOUNDARY 


0.75k 


1.00k 


1.18k 


1.30k 


- FINITE  ELEMENT  SOLUTION 

_ _ 1 _ I _ I _ L 


NOTE:  THE  APEX  OF  THE  CURVES  REPRESENTS  THE  BOUNDARY  BETWEEN 
THE  PLASTIC  AND  ELASTIC  REGIONS.  IN  COMPARING  CURVES 
IT  IS  NECESSARY  TO  EXAMINE  CURVES  WITH  THE  SAME  p/a. 


CLOSED  FORM  SOLUTION 
CPRAGER  AND  HODGE,  1951) 
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INCREMENT  II 
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- —CLOSED  FORM  SOLUTION 

(PRAGER  S  HODGE,  1951) 

- FINITE  ELEMENT  SOLUTION 
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•  -  INCREMENT  1,  p  =  0.75k 
O-  INCREMENT  2,  ITERATION  3,  p  =  1.00k 
A-  INCREMENT  3,  ITERATION  6,  p  =  1.18k 
□-  INCREMENT  4,  ITERATION  6,  p  =  1.30k 
O-  INCREMENT  5,  ITERATION  6,  p  =  1.39k 
p  =  RADIUS  TO  ELASTO-PLASTIC  BOUNDARY 


FIG.  12-  DISTRIBUTION  OF  RADIAL  STRESS  FOR  A 
THICK-WALLED  CIRCULAR  TUBE 
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FIG.  13  -  DISTRIBUTION  OF  AXIAL  STRESS  FOR  A 
THICK-WALLED  CIRCULAR  TUBE 
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p/a 


FIG.  14  -  RADIAL  DISPLACEMENT  U  (a) ,  U  (b)  VS.  RADIUS  p  OF  ELASTIC- 
PLASTIC  BOUNDARY  FOR  A  THICK -WALLED  CIRCULAR  TUBE 
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circular  opening  shown  in  Fig.  15  was  analyzed.  The  definition 
of  the  problem  and  the  finite  element  idealization  are  shown  in 
Fig.  15.  The  analysis  was  conducted  by  applying  pressures  on  the 
outer  boundary.*  The  boundary  pressures  which  have  vertical  and 
horizontal  components  were  applied  in  14  increments.  The  percen¬ 
tage  of  the  total  boundary  pressure  applied  for  each  increment  is 
tabulated  in  Fig.  15. 

It  was  experienced  that  the  solution  convergence  was  slow  as 
compared  with  the  case  for  the  thick-walled  circular  tube  des¬ 
cribed  above.  This  convergence  problem  was  also  noted  by  Baker, 
et  al.  (1969).  For  the  purpose  of  speeding  convergence,  the 
following  over-relaxation  technique  was  adopted. 

For  every  iteration,  the  excess  stresses  (Ao"}  as  computed  by 
Eq.  (22),  which  are  to  be  balanced  by  a  set  of  body  forces,  are 
multiplied  by  an  over-relaxation  factor  to  bring  the  element  to 
the  point  below  the  yield  surface.  The  over-relaxation  factor, 

Ro  (p) ,  which  will  ensure  convergence  of  the  solution,  was  found 
to  be  in  the  range  between  1.0  and  1.5.  Its  value  was  selected 
in  accordance  with  the  following  procedure. 


*  In  the  semi-annual  technical  report,  this  problem  was 

analyzed  with  different  loading  conditions.  A  comparison 
of  the  results  is  presented  in  the  Appendix  D. 
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FIG.  15  -  FINITE  ELEMENT  MESH  FOR  AN  ELASTO- PLASTIC  ANALYSIS  OF 
A  CIRCULAR  OPENING 
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(1)  For  every  iteration  f  ,  which  is  an  index  showing  where 
the  element  was  situated  outside  the  yield  surface,  is 
calculated  according  to 

fr  =  [£(o')  -  k]  /  (k  -  aJx)  (27) 

(2)  Define  the  current  state  of  stress  by  the  following 
equation 

{o}  *  {o')  -  (A0")1  •  p  (28; 

p  is  to  be  determined.  Express  [f(o)  -  k]  in  terms  of  p. 

(3)  (a)  If  f  >_  15%,  determine  p  on  the  basis  of  the  following 

equation: 

[f (0)  -  k]  -  [k  -  f(o') J  (29) 

(b)  If  15%  <  ff  ^  10%,  a  value  of  p  was  selected  such  that 

f  (0)  -  k  <  0.75  [k  -  f(0')]  (30) 

(c)  If  f  <  10%,  a  value  of  p  was  selected  such  that 

f(0)  -  k  <  0.5  [k  -  f (0 ' ) J  (31) 

Using  this  technique  for  increasing  the  rate  of  convergence, 
six  to  eight  iterations  were  still  needed  for  a  loading  incre¬ 
ment.  The  results  of  the  analysis  together  with  those  pre¬ 
sented  by  Reyes  (1966)  and  Baker,  et  al  (1969)  are  shown  in 
Figures  16  and  17  .  Figure  16  shows  vertical  (tangential)  and 
horizontal  (radial)  stresses  along  a  horizontal  section  through 
the  opening.  All  three  methods  of  analysis  give  the  same 
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FIG.  16  -  VERTICAL  AND  HORIZONTAL  STRESSES  ALONG  HORIZONTAL 


PRESENT  STUDY 

REYES  0  966) 

INFERRED  FROM  BAKER, 
SANDHU,  SHIEH  0  969) 
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FIG.  17  -  DEFORMATION  ALONG  CAVITY  FACE  OF  A  CIRCULAR  OPENING  AS 
COMPUTED  BY  ELASTIC  AND  ELASTIC -PLASTIC  ANALYSIS 
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distribution  of  the  horizontal  stress.  The  tangential  stress 
distribution  indicates  that  the  three  methods  of  analysis  predict 
approximately  the  same  extent  of  the  yield  zone.  The  peak  tan¬ 
gential  stress,  however,  differs  by  approximately  251.  This 
difference  is  thought  to  be  due  to  two  causes:  (i)  differences 
in  the  degree  of  convergence  obtained,  and  (ii)  different  methods 
of  computing  the  axial  stress.  As  discussed  earlier,  the  yield 
criterion  is  assumed  to  be  a  function  of  the  sum  of  three  normal 
stresses  and  the  second  invariant  of  the  stress  deviation.  There¬ 
fore,  different  values  of  the  axial  stress  result  in  different 
yield  surfaces.  The  differences  in  the  methods  utilized  by  Reyes 
(1966),  Baker,  et  al .  (1969)  and  that  utilized  in  this  study  to 
calculate  the  axial  stress  are  briefly  discussed  in  the  following 
paragraphs . 

Baker,  et  al.  (1969)  used  the  following  equation  to  compute  the 

axial  stress ,  a, . 

z 

In  the  elastic  range,  the  axial  stress  is  given  by 

oz  =  v  (CTj  +  o3)  (32) 

where  and  o3  are  major  and  minor  principal  stress,  respectively, 
in  the  x-y  plane,  and  v  =  Poisson's  ratio. 

In  the  plastic  range,  the  axial  stress  is  expressed  by 

°z  =  I  +  a3)  '  J  Ccr1  -  a3)  Sin  <f>  (33) 


WOODWARD-LUNDGREN  S.  ASSOCIATES 


-50- 


where  4»  =  the  angle  of  internal  friction. 


Comparison  of  equation  (32)  and  (33)  shows  that  the  values  of 

a  computed  by  the  two  equations  will  not  be  the  same  for  the 
z 

cases  in  which  v  f  j  and  <}>  f  0  indicating  that  there  will  be 

an  abrupt  change  on  the  value  of  a  as  the  element  undergoes 

from  the  elastic  state  to  the  plastic  state.  The  total  stress 

is  used  at  each  computation.  A  different  method  was  used  by 

Reyes  (1966)  to  calculate  a  .  For  every  increment  of  loading, 

z 

o7  was  first  calculated  using  the  equation  similar  to  equation 
(16)  assuming  that  the  element  was  in  the  elastic  range. 

(r)  .  „  (r-1)  +  v  (A  (r)  +  a  (r),  (34) 

z  z  j\.  y 

If  the  yield  function  was  exceeded  and  the  element  was  in  yield 
in  the  previous  load  increment,  the  new  state  of  stress  was  com¬ 
puted  by 


(r)  = 


(r-1) 


ij 


ij 


+  2G(Ae ■ 


ij 


a6 


-  bo  ■ 


(r) 


(35) 


where  a—  are  stress  components,  G' is  shear  modulus,  Ae  —  are 
the  computed  strain  increments,  6^  is  the  Kronecker  delta,  and 
a,  b  are  constants  to  be  determined,  Equation  35  represents  a 
set  of  equations  for  each  stress  component.  These  equations 
were  solved  for  the  constants  a  and  b,  and  the  new  state  of 
stress  ^  by  successive  substitutions. 


The  present  study  used  a  different  approach  to  evaluate  a  . 

The  loading  is  applied  incrementally,  and  for  each  load  increment 
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the  change  in  axial  stress  Acz  is  calculated  in  accordance  with 
equations  16  or  17  depending  upon  whether  the  element  is  in  the 
elastic  or  plastic  range. 

The  vertical  and  horizontal  components  of  the  deflection  at  the 
face  of  the  opening  are  presented  in  Fig.  17.  It  may  be  seen 
that  the  elasto -plastic  solution  indicates  an  increase  in  both 
vertical  and  horizontal  deflections  toward  the  opening  as  com¬ 
pared  with  the  elastic  solution.  The  horizontal  deflection  for 
the  solution  obtained  by  Reyes  (1966)  does  not  appear  to  be  reli¬ 
able.  The  comparison  of  the  vertical  deflections  indicates  ex¬ 
cellent  agreement  between  the  results  obtained  by  Reyes  (1966) 
and  the  present  study. 

III.  Analysis  of  a  Rectangular  Underground  Opening 
To  demonstrate  the  usage  of  the  combined  computer  program  devel¬ 
oped  in  which  all  three  types  of  analysis  -  no  tension,  joint 
perturbation  and  elasto-plastic  analysis  are  considered,  a  hypo¬ 
thetical  rectangular  opening  located  at  the  depth  of  1000  ft  was 
analyzed.  The  finite  element  idealization  is  shown  in  Fig.  18. 

A  horizontal  joint  is  assumed  to  exist  50  ft  above  the  crown  of 
the  opening.  Initial  stresses  were  assumed  to  increase  with 
depth  in  accordance  with  the  gravity  stress  field.  The  vertical 
stress  is  assumed  equal  to  depth  times  unit  weight  of  the  rock 
(144  pcf  assumed)  and  the  horizontal  stress  equal  to  v/(l-v) 
times  the  vertical  stress.  The  problem  is  defined  in  Fig.  18. 

The  excavation  for  the  opening  was  simulated  by  applying  the 
boundary  pressures  in  four  increments  as  shown  in  Fig.  18. 
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FIG.  18  -  FINITE  ELEMENT  MESH  FOR  ANALYSIS  OF  A  RECTANGULAR  OPENING 
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The  results  of  the  analysis  are  presented  in  Figs.  19  through 
25.  Figs.  19  and  20  show,  respectively,  the  distribution  of 
the  normal  and  tangential  stress  along  the  horizontal  joint 
for  both  elastic  and  combined  solution.  Some  readjustment  of 
stresses  may  be  noted  near  the  centerline  of  the  opening  re¬ 
sulting  from  the  yielding  of  the  joint  elements.  It  should, 
however,  be  recognized  that  only  a  few  joint  elements  have 
yielded  and  there  are  no  large  movements  along  the  joint  since 
the  joint  is  intact  over  the  majority  of  its  length.  Figs.  21 
and  22  show,  respectively,  the  distribution  of  the  major  prin¬ 
cipal  stress  in  the  vicinity  of  the  opening  obtained  from  the 
elastic  and  combined  solution;  Figs  23  and  24  are  for  the  dis¬ 
tribution  of  the  minor  principal  stress  obtained  from  the  elastic 
and  combined  solution.  It  may  be  noted,  by  comparing  with  the 
elastic  solution,  that  some  readjustment  of  stresses  occurs  in 
the  vicinity  of  the  opening.  The  magnitude  and  region  of  tensile 
stresses  have  been  reduced  by  the  stress  transfer  process  to  simu¬ 
late  the  no  tension  rock  characteristic.  Fig.  25  indicates  the 
development  of  the  plastic  zone  in  the  vicinity  of  the  opening. 

Summary  of  Achievement  for  Phase  lc 

The  combined  computer  program  developed  in  Phase  lb  was  used  to 
analyze  certain  illustrative  problems.  Satisfactory  agreement 
was  obtained  when  the  results  of  these  analyses  were  compared 
with  closed  form  solutions  and  the  solutions  obtained  by  other 
researchers.  Using  this  combined  computer  program,  complex  geo¬ 
metries  and  heterogeneous  rock  conditions  and  properties  can  be 
accounted  for. 
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FIG.  20  -  DISTRIBUTION  OF  TANGENTIAL  STRESS  ALONG  HORIZONTAL  JOINT 
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FIG.  22  -  DISTRIBUTION  OF  MAJOR  PRINCIPAL  STRESS  AROUND 
A  RECTANGULAR  OPENING  (COMBINED  FINITE  ELEMENT 
SOLUTION) 
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FIG.  24  -  DISTRIBUTION  OF  MINOR  PRINCIPAL  STRESS  AROUND  A 
RECTANGULAR  OPENING  (COMBINED  FINITE  ELEMENT 
SOLUTION) 
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DEVELOPMENT  OF  PLASTIC  ZONES  AROUND  A  RECTANGULAR 
OPENING  (COMBINED  FINITE  ELEMENT  SOLUTION) 
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ANALYSES  OF  LABORATORY  MODELS  AND  CASE  HISTORIES  (Phase  7) 

The  Phase  2  work  consists  of  the  analyses  of  a  laboratory  model 
study  and  a  case  history  study  using  the  combined  computer  program 
developed  in  Phase  1.  The  results  of  these  analyses  were  compared, 
with  observed  performance.  Selection  of  case  histories  and  labora¬ 
tory  models  w as  made  on  the  basis  of  availability  of  adequate 
information  on  material  properties ,  geological  conditions  and 
performance  to  satisfactorily  model  the  problem. 

The  possibility  of  utilizing  the  results  of  tests  on  well- 
instrumented  laboratory  models  of  unlined  openings  in  rock-like 
materials  e.g.  ,  Heuer  and  Hendron  (1971)  ,  Haas  and  Clark  (1970)  , 
and  Rosenblad  (1971)  for  evaluating  the  ability  of  the  program 
developed  in  Phase  1  to  predict  performance  was  also  studied. 

The  geologic  conditions  and  performance  of  the  Morrow  Point 
Underground  Powerplant  and  Straight  Creek  Tunnel  Pilot  Bore  were 
studied  to  determine  their  suitability  for  analysis. 

Preliminary  studies  indicated  that  the  Morrow  Point  Underground 
Powerplant  and  model  tests  reported  by  Heuer  and  Hendron  (1971) 
were  suitable  for  evaluating  the  ability  of  the  computational 
techniques  developed  in  Phase  1  to  predict  stresses,  strains, 
and  deformations  in  a  rock  mass  surrounding  an  excavation. 


Analysis  of  Laboratory  Model 
Description  of  Model  Study 

Heuer  and  Hendron  (1971)  conducted  a  series  of  model  tests  on 
unlined  openings  in  a  rock-like  material.  The  geometry  of  the 
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model  block  tested  is  shown  in  Fig.  26.  The  model  was  constructed 
to  represent  a  segment  of  a  long  tunnel.  The  test  was  conducted 
under  plane  strain  conditions  which  were  achieved  by  adjusting 
the  axial  stress  in  the  longitudinal  direction  during  the  test 
to  null  any  longitudinal  displacement.  The  model  block  consisted 
of  two  halves  assembled  as  shown  in  Fig.  26.  The  behavior  of  the 
model  was  monitored  by  internal  foil  strain  gage  rosettes  installed 
on  the  block  interface  at  various  radii.  Circumferential  strain 
gages  were  also  installed  on  the  tunnel  wall.  Two  sets  of  dia¬ 
metrical  extensometers  were  also  used  to  measure  the  closure  of 
the  opening.  The  loads  applied  to  the  model  consisted  of  hori¬ 
zontal  and  vertical  stresses  applied  at  the  boundary  as  shown  in 
Fig.  26.  Tests  were  conducted  under  three  ratios  of  horizontal 
to  vertical  stress,  (a)1/av)  »  l-0»  -75  and  .25.  Only  the  tests 
with  <Jj1/ov  =  *25  are  analyzed  in  this  study. 

The  material  used  for  construction  of  the  model  blocks  was  a 
water/plaster/sand  mixture  in  the  ratio  of  1. 2/1/9.  The  sand 
was  a  fine  Sangamon  sand.  A  summary  of  the  test  results  obtained 
from  triaxial  tests  on  the  material  is  presented  in  Figs.  27  to 
29.  From  Fig.  27,  the  material  constants  used  in  Equation  (5) 
for  the  yield  function  were  calculated  to  be  a  =  0.262  and  k  =  139  psi 
respectively.  The  tensile  strength  of  the  material  was  determined 
to  be  35  psi  as  compared  with  the  unconfined  compressive  strength 
of  600  psi.  Fig.  28  indicates  the  nonlinear  nature  of  the  axial 
stress  -  strain  curves  at  different  confining  pressures.  The  stress- 
strain  curves  show  that  the  linear  approximation  is  valid  only  up 


WOODWARD-LUNDGREN  &  ASSOCIATES 


-63- 


14" 


FIG. 


26  -  GEOMETRY  OF  MODEL  BLOCK 

(AFTER  HEUER  AND  HENDRON ,  1971) 


WOODWARD-IUNDGREN  &  ASSOCIATES 


29.6°  qu  =  600  PSI 

133  TENSILE  STRENGTH 

0.262  k  =  189  PSI  _ 
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FIG.  27  -  SUMMARY  OF  TRIAXIAL  COMPRESSION  TEST  ON  ROCK-LIKE  MATERIAL 
(AFTER  HEUER  AND  HENDRON ,  1971) 


o  r2  .4  .6  , 8  1.0 

AXIAL  STRAIN,  6  C%) 


FIG.  28  -  AVERAGE  STRESS -STRAIN  CURVES  AT  DIFFERENT  CONFINING 
PRESSURES  -  TRIAXIAL  COMPRESSION  SPECIMENS  (AFTER 
HEUER  AND  HENDRON,  1971) 
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to  the  principal  stress  difference  (a  -  a3)  of  600  to  800  psi 
depending  on  the  confining  pressure.  Fig.  29  shows  the  variation 
of  the  Poisson's  ratio  with  the  principal  stress  difference.  It 
may  be  noted  that  the  Poisson's  ratio  increases  with  an  increase 
in  the  principal  stress  difference,  the  increase  depending  on 
the  confining  pressure.  This  behavior  indicates  that  the  mate¬ 
rial  dilates  during  shearing.  The  value  of  the  Poisson's  ratio 
as  measured  from  the  test  varied  from  0.14  to  0.67.  Values  above 
0.5  indicate  that  the  material  is  not  an  ideal  linear  isotropic 
elastic  material. 

Behavior  of  the  Model 

The  behavior  of  the  model  during  loading  for  the  case  analyzed 
(ah/av  -  1/4)  is  described  in  a  later  section  in  which  the  com¬ 
puted  and  observed  behavior  is  compared.  Only  the  development 
of  fractures  is  discussed  here.  Fig.  30  shows  the  pattern  of 
fractures  developed  at  ay  =  800  psi.  During  testing,  the  frac¬ 
tured  wedges  at  the  springlines  shown  in  Fig.  30  became  visible 
when  the  applied  load  was  incremented  from  ay  =  690  psi  to 
ay  =  765  psi.  As  the  stress  level  was  increased,  more  fractures 
formed  and  the  wedges  moved  into  the  opening.  Another  signifi¬ 
cant  set  of  fractures  is  seen  in  Fig.  30  extending  away  from  the 
opening  back  into  the  model.  These  fractures  appear  to  be  exten¬ 
sions  of  the  fractures  which  formed  the  wedges  at  the  springline. 
At  the  crown  and  invert,  high  tensile  circumferential  strains  were 
measured  at  applied  stress  levels  above  ctv  =  400  psi.  Heuer  and 
Hendron .  (1971)  hypothesized  that  tension  cracks  might  have  formed 
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FIG.  29  -  VARIATION  OF  POISSON'S  RATIO  WITH  PRINCIPAL  STRESS  DIFFERENCE 
TRIAXIAL  COMPRESSION  TESTS  (AFTER  HEUER  AND  HENDRON,  1971) 
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FIG.  30  -  FRACTURES  DEVELOPED  IN  MODEL  TEST 
(oh/av  =  1/4)  (AFTER  HEUER  AND 
HENDRON,  1971) 
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during  testing.  However,  no  cracks  were  visible  after  testing, 
It  was  assumed  that  the  cracks  were  open  during  the  test  and 
then  closed  upon  removal  of  load. 

Analysis  of  Model  Study 

Idealization  of  the  Model  -  Because  of  the  symmetry  of  the 
model  with  respect  to  its  axes,  as  shown  in  Fig.  26,  it  is 
only  necessary  to  analyze  a  quadrant  of  the  cross-section. 

Fig.  31  shews  the  finite  element  idealization  of  the  model. 

Material  Properties  Used  in  the  Analysis  -  As  described  previ¬ 
ously,  the  material  tested  shows  some  nonlinear  stress-strain 
behavior  with  respect  to  the  elastic  modulus  and  the  Poisson's 
ratio.  For  the  purpose  of  analysis,  it  was  necessary  to  use  a 
single  value  of  the  modulus  and  the  Poisson's  ratio.  The  value 
of  the  modulus  was  determined  from  Fig.  36  to  be  833,000  psi. 
The  Poisson's  ratio  of  0.14  was  used  in  the  analysis. 

The  incremental  analysis  was  performed  to  simulate  the  actual 
loading  path  =  1/4).  The  incremental  loads  applied  are 

shown  in  Fig.  31.  Two  Cases  were  analyzed  under  the  assumption 
that  the  material  behaved  as  an  elasto -plastic  material  and  had 
a  limited  tensile  strength  of  35  psi.  The  first  case  (Case  1) 
assumed  that  there  was  no  fractured  wedge  formed  during  testing 
The  second  case  (Case  2)  assumed  that  a  fractured  wedge  formed 
at  the  springline  during  testing  ,-utd  separated  from  the  main 
body  of  the  model.  The  details  of  these  cases  are  described 
in  the  following  sections. 
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Presentation  and  Discussion  of  Results 

The  results  of  analysis  together  with  some  of  the  measured 
behavior  are  shown  in  Figs.  32  to  38.  Figs.  32  and  33  illustrate 
the  development  of  plastic  and  tensile  regions  at  three  levels 
of  the  axial  stress  400,  600,  and  800  psi.  The  observed  frac¬ 
tures  are  also  shown  in  Figs.  32  and  33.  On  comparing  Figs.  30 
and  32,  it  may  be  noted  that  the  development  of  the  plastic 
region  in  Fig.  32  is  similar  to  the  shear- fractured  wedges 
developed  at  the  springline  as  shown  in  Fig.  30.  At  oy  =  800  psi, 
a  large  plastic  region  developed  back  into  the  model,  Fig.  32, 

The  location  of  the  actual  shear  fracture,  also  shown  in  Fig.  32, 
is  along  the  edge  of  the  plastic  region.  It  was  observed  in  the 
model  test  that  a.  region  of  tensile  stress  developed  at  the  crown 
and  invert.  However,  there  was  little  failure  in  this  zone  indi¬ 
cating  that  the  tensile  stresses  were  below  the  tensile  strength 
of  the  material.  An  elastic  analysis  indicates  that  a  large 
tensile  region  developed  with  a  tensile  stress  higher  than  the 
tensile  strength  at  the  crown  and  invert  (Fig.  34).  An  elasto- 
plastic  analysis  indicates  that  a  much  smaller  zone  of  tension 
cracks  would  develop  if  the  material  exhibits  elasto -plastic 
behavior.  The  results  of  the  elasto-plastic  analysis  appears 
to  agree  with  observed  behavior. 

In  order  to  simulate  the  propagation  of  fractures  and  their 
effects  on  the  behavior  of  the  model,  an  analysis  was  performed 
(Case  2)  after  removing  the  wedge  that  formed  under  a  vertical 
stress  of  600  psi  as  indicated  by  the  results  of  the  analysis 
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FIG.  33  -  RESULTS  OF  ANALYSIS  OF  MODEL  TEST,  ah/ov  =  1/4,  SHOWING  DEVELOPMENT  OF 
PLASTIC  REGIONS,  CASE  2 
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RESULTS  OF  ANALYSIS  OF  MODEL  TESTS 
SHOWING  TENSILE  REGION,  CASE  2 
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DISTRIBUTION  OF  HORIZONTAL  STRAINS  ALONG  6  =  86.25 
ANALYSIS  OF  MODEL  TESTS 
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shown  in  Fig.  32.  Under  these  conditions,  the  plastic  zone 
extended  behind  the  wedge  and  propagated  into  the  model  block 
as  shown  in  Fig.  33.  The  direction  in  which  the  plastic  zone 
developed  as  the  loads  increased  was  similar  to  the  observed 
propagation  of  fractures.  Fig.  33  shows  the  extend  of  the 
plastic  zone  as  the  loads  were  increased  to  the  vertical  stress 
of  800  psi.  It  may  be  noted  that  the  observed  shear  fracture 
lies  in  the  middle  of  the  plastic  zone  indicating  that  the 
elastic- analys is  can  provide  a  good  indication  with  respect  to 
the  development  of  the  critical  zone  in  the  vicinity  of  the 
opening . 

The  tensile  stress  region  in  the  final  stress  distribution  and 
the  zone  where  the  major  principal  tensile  stress  was  greater 
than  the  tensile  strength  during  the  entire  loading  process  is 
shown  in  Fig.  35.  Comparisons  of  tensile  regions  and  the  regions 
with  values  of  tensile  stress  greater  than  tensile  strength  during 
loading  for  Cases  1,  2  and  the  elastic  analysis  (Figs.  32,  33  and  34, 
respectively),  indicate  the  following  features:  (i)  The  Case  2 
analysis  shows  that  a  larger  region  with  a  tensile  stress  greater 
than  the  tensile  strength  (35  psi)  developed  during  loading. 

(ii)  The  "no  tension"  analysis  indicates  that  upon  subsequent 
loading  the  tensile  region  near  the  crown  and  invert  was  recom¬ 
pressed.  (iii)  The  tensile  region  was  shifted  from  the  crown  and 
invert  to  the  region  closer  to  the  fractures  at  the  springline 
indicating  that  the  development  of  fractures  at  the  springline 
caused  a  redistribution  of  stresses  in  the  rock  surrounding  the 
opening . 
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In  order  to  compare  with  the  strains  measured  in  the  model  tests, 
element  strains  were  computed  from  the  nodal  point  displacements. 
These  strains  were  assumed  to  occur  at  the  centroid  of  the  ele¬ 
ment.  The  computed  vertical  and  horizontal  strains  along  the 
line  (0  =  86.25°)  close  to  the  horizontal  axis  (0  =  90°)  for 
Cases  1  and  2  are  shown  in  Figs.  36  and  37,  respectively,  for 
the  stress  levels  between  o =  600  psi  and  ay  =  800  psi.  For 
the  purpose  of  comparison,  the  measured  tangential  and  radial 
strains  along  the  horizontal  axis  are  also  plotted  in  Figs.  36 
and  37.  Because  of  the  relatively  large  scatter  in  the  measure¬ 
ments  obtained  from  the  electrical  strain  gages,  Heuer  and  Hendron 
(1971)  concluded  that  the  strain  measurements  would  only  give  a 
qualitative  indication  of  the  response  of  the  model.  Comparison 
of  the  measured  and  calculated  behavior  shows  that  at  the  low 
stress  levels,  they  are  in  good  agreement.  Figs.  36  and  37. 
However,  as  the  stress  level  increased  and  the  fractured  wedges 
formed,  the  strains  increased  to  a  level  much  greater  than  that 
predicted  from  the  Case  1  analysis  in  which  the  model  block  was 
assumed  to  act  as  continuum.  The  results  of  the  Case  2  analysis 
indicated  a  strain  pattern  and  magnitude  in  approximate  agreement 
with  measured  values.  This  indicates  that  the  formation  of  frac¬ 
tured  wedges  at  the  springline  has  a  great  effect  on  the  behavior 
of  the  model  block.  Discrepancies  between  the  measured  strains 
and  the  strains  computed  by  Case  2  analysis  at  the  higher  stress 
levels  may  in  large  part  be  attributed  to  the  formation  of  the 
second  set  of  fractures  extending  back  into  the  model  block  away 
from  the  opening. 
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Fig.  38  shows  the  measured  and  computed  radial  displacements 
at  a  point  on  the  horizontal  axis  (0  =  90°,  r/a  =  1.65).  Both 
the  measured  and  computed  radial  displacements  indicate  the 
point  moved  outward  from  the  opening  as  the  loads  increased. 

The  measured  displacements  are  about  three  times  larger  than 
those  computed  from  the  elastic  or  elasto-plastic  analysis 
from  the  low  stress  level  up  to  =  600  psi.  The  large  in¬ 
crease  .in  displacements  for  a  '  £600  psi  may  be  attributed  to 
the  development  of  the  fractured  wedges  and  the  fractures  back 
into  the  model  block.  The  Case  2  analysis  showed  a  similar 
increase  in  radial  displacement  due  to  the  formation  of  the 
fractured  wedges  at  the  springline.  Discrepancies  between  the 
measured  and  calculated  displacements  at  higher  stress  levels 
may  also  be  attributed  to  the  formation  of  the  second  set  of 
fractures  back  into  the  model  block. 

Analysis  of  the  model  tests  indicates  that  the  elasto-plastic 
analysis  provides  valuable  information  for  evaluating  the  criti¬ 
cal  zone  in  an  intact  rock  surrounding  an  opening. 

After  significant  fractures  develop  in  the  rock  surrounding  an 
opening,  the  rock  becomes  discontinuous  along  the  fractures  and 
strain  displacement  distributions  within  the  rock  could  not  be 
accurately  determined  by  the  elasto-plastic  analysis.  This  is 
due  to  the  inability  of  the  analysis  to  simulate  (1)  the  prop¬ 
agation  of  shear  fractures  in  the  rock,  and  (2)  the  nonlinear 
stress-strain  behavior  and  dilatancy  of  the  rock-like  material 
during  shearing. 
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Although  the  analysis  cannot  simulate  accurately  the  behavior 
of  the  model  after  the  formation  of  shear  fractures  surrounding 
the  opening,  some  indication  can  be  obtained  with  regard  to 
propagation  of  fractured  zones. 

Morrow  Point  Underground  Powerplant  Excavation 
Description 

The  Morrow  Point  Underground  Powerplant  was  constructed  by  the 
U.S.  Bureau  of  Reclamation  on  the  Gunnison  River  some  20  miles 
east  of  Montrose,  Colorado  (Dodd,  1967  and  Brown,  et  al.  1971). 

The  powerplant  chamber  is  206  ft  long  and  57  ft  wide  with  a 
height  ranging  from  65  to  134  ft  and  is  located  about  400  ft 
below  the  ground  surface.  The  plan  of  the  powerplant  and  other 
adjacent  structures  is  shown  in  Fig.  39.  A  cross-section  of  the 
powerplant  chamber  along  A-A’ (line  4  +  12  ft)  is  shown  in  Fig.  40. 
It  may  be  noted,  from  Figs.  39  and  40,  that  the  powerplant  chamber 
is  situated  behind  a  steep  valley  rock  wall.  The  crown  of  the 
chamber  on  the  river  side  (b  line  wall,  Fig.  40)  is  about  200  ft 
behind  the  ground  surface. 

The  powerplant  is  located  entirely  within  the  Pre-Cambrian 
metamorphic  rocks.  The  composition  of  the  metamorphic  rock 
consists  predominantly  of  mica  schist  and  quartz -mica  schist 
with  some  biotite  schist  and  pegmatite.  The  bedding  strikes 
nearly  normal  to  the  powerplant  alignment  and  dips  upstream 
at  angles  varying  from  15°  to  60°,  averaging  about  35°.  There 
are  two  distinct  shear  zones  intersecting  the  powerplant  area. 

The  lower  zone  (shear  zone  A)  strikes  N  40°  W,  dips  32°  E,  and 
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FIG.  39  -  PLAN  AND  LOCATION  OF  MORROW  POINT  POWERPLANT  (AFTER  DODD,  1967) 
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FIG.  40  -  CROSS-SECTION  A-A’  (LINE  4  +  12  FT)  OF  THE  POWERPLANT 
CHAMBER,  MORROW  POINT  POWERPLANT 


WOODWARD-IUNDGREN  &  ASSOCIATES 


-85- 


consists  of  a  one  to  five  ft  zone  of  fau  t  gouge  and:  fractured 
rock.  The  upper  zone  (shear  zone  B)  has  a  strike  of  N  20°  E, 
a  dip  of  23°  E  and  consists  of  a  one  to  three  ft  zone  of  fault 
gouge  and  fractured  rock.  The  orientations  of  the  three  major 
joint  sets  that  intersect  the  chambers  are:  strike  N  63°  W  and 
dip  82°  SW ;  strike  N  36°  E,  and  dip  80°  W;  strike  N-S  and  dip 
43°  E. 


Five  basic  rock  types  are  present  in  the  chamber.  Their  per¬ 
centages  as  exposed  at  the  rock  surface  of  the  chamber  are 
shown  below: 


Type 

I 

Biotite  Schist 

15 

Type 

II 

Mica  Schist 

20 

Type 

III 

Quartz -Mica  Schist  or  Micaceous 
Quartzite 

40 

Type 

IV 

Quartzite 

20 

Type 

V 

Pegmatite 

5 

Significant  Features  of  the  Behavior  of  the  Powerplant  Excavation 
The  crown  of  the  chamber  was  first-  excavated  during  a  three  month 
period  from  May  1965  through  July  1965.  The  excavation  of  the 
remaining  chamber  continued  from  August  1965  through  March  1966. 
When  the  excavation  reached  the  El.  6793  bench,  48  ft  below  the 
crown,  some  initial  movement  of  the  a-line  rock  wall  occurred. 
When  the  powerplant  was  excavated  to  El.  6748,  92  ft  below  the 
crown,  an  accelerated  inward  movement  of  1.5  inches  was  measured 
at  line  4+12  ft.  As  the  excavation  proceeded,  the  a-line  rock 
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wall  continuously  moved  inward  but  at  a  diminishing  rate. 
Significant  rock  movement  ceased  by  the  end  of  March  1966 
(Brown,  Morgan  and  Dodd,  1971).  The  maximum  observed  inward 
movement  of  the  a-line  wall  at  line  4  +  12  ft  was  2.5  inches. 

Fig.  41  shows  the  history  of  the  observed  rock  movements  at 
line  4  +  12  ft  reported  by  Brown,  et  al.  (1971). 

Brown  et  al .  (1971)  hypothesized  that  the  anomalous  rock 
behavior  along  the  a-line  rock  wall  might  involve  a  mass 
movement  of  a  rock  wedge  shown  in  Figs.  42  and  43.  It  was 
hypothesized  that  shear  zones  A  and  B  intersect  at  an  average 
distance  of  130  ft  behind  the  a-line  wall,  forming  a  rock 
wedge.  In  addition,  two  planes  of  failure  must  exist  for  the 
wedge  to  move  into  the  opening.  One  is  the  upper  shear  failure 
plane  with  an  apparent  dip  of  17°  intersecting  the  chamber  near 
the  springline  of  the  rock  arch.  The  other  is  essentially  ver¬ 
tical  and  intersects  the  chamber  between  line  2  +  15.5  ft  and 
the  east  end  wall  of  the  chamber.  A  cross-section  through  the 
chamber  at  line  4  +  12  ft  showing  the  rock  wedge  is  illustrated 
in  Fig .  43 . 

To  stabilize  rock  movements  along  the  a-line  wall,  the  following 
four  steps  were  taken  (Brown,  et  al . ,  1971).  (1)  Installation 

of  nine  1.69  inch  diameter  deformed  reinforcement  anchor  bars 
in  drill  holes  crossing  shear  zone  A  from  the  hanging  wall  to 
the  foot  wall.  (2)  Installation  of  27  additional  long  rock  bolts 
in  the  a-line  wall  between  the  1-line  and  3-line.  The  bolts  were 
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FIG.  43  -  SCHEMATIC  DIAGRAM  OF  A-LINE  ROCK  WEDGE  MOVEMENT 
(AFTER  BROWN,  ET  AL.  1971) 
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placed  diagonally  in  plan  across  the  vertical  zone  of  incipient 
failure  into  sound  rock.  (3)  Installation  of  25  post-tensioned 
tendons  from  the  a- line  rock  wall  through  shear  zone  A  to  the 
drainage  tunnel.  (4)  Installation  of  10  pairs  of  flat  jacks  in 
the  first  stage  concrete  opposite  the  exposed  face  of  the  rock 
wedge  to  provide  a  restraining  force  between  the  concrete  struc¬ 
ture  and  rock. 

Idealization  of  the  Powerplant  Excavation 

a.  Finite  Element  Idealization  -  As  discussed  previously,  the 
observed  maximum  rock  movements  occurred  at  line  4  +  12  ft. 
Therefore,  a  cross-section  through  the  opening  at  this  location 
(Fig.  43)  was  chosen  for  analysis.  A  finite  element  idealization 
of  the  section  is  shown  in  Fig.  44.  The  essential  features  in 
the  idealization  are  the  steep  sloping  valley  ground  surface, 
shear  zones  A  and  B  and  two  other  incipient  failure  planes.  The 
steep  sloping  ground  surface  makes  estimating  the  initial  in-situ 
stresses  difficult.  As  described  previously,  both  the  rock  wedge 
and  the  movements  associated  with  it  are  three-dimensional  in 
nature.  However,  for  the  purposes  of  analysis,  it  was  necessary 
to  idealize  it  as  a  two-dimensional  plane  strain  problem.  Shear 
zones  A,  B  and  the  other  two  incipient  failure  planes  are  shown 
in  Fig.  43  as  distinct  discontinuities  extending  in  the  direction 
perpendicular  to  the  cross-section.  These  discontinuities  are 
idealized  with  Goodman's  one -dimensional  joints. 

b.  Material  Properties  -  Extensive  testing  programs  were  con¬ 
ducted  by  the  Bureau  of  Reclamation  (1965)  to  determine  rock 
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properties  at  the  Morrow  Point  Dam  site.  These  included  direct 
shearing  and  sliding  friction  both  in  the  field  and  the  laboratory, 
triaxial  tests  and  field  jacking  tests  of  foundation  rock.  The 
strength  and  elastic  properties  of  the  rock  above  and  below  shear 
zone  A  used  in  the  analysis  are  summarized  in  Table  1.  As  des¬ 
cribed  previously,  both  shear  zones  A  and  B  were  idealized  by 
one-dimensional  joint  elements.  The  properties  of  these  shear 
zones  are  approximated  by  the  normal  and  shear  joint  stiffnesses 
(Goodman,  1969)  which  are  functions  of  normal  and  tangential 
deformability  of  the  shear  zones  and  their  geometries.  No  data 
were  available  for  the  normal  and  shear  stiffnesses  of  the  shear 
zones.  Therefore,  two  sets  of  joint  stiffness  were  utilized  to 
parametrically  study  the  influence  of  the  deformability  of  the 
sheai  zones  on  the  behavior  of  the  powerplant  excavation.  The 
strength  parameters  assumed  for  the  shear  zones  are  also  presented 
in  Table  1. 

c.  Initial  State  of  Stress  at  the  Powerplant  -  The  initial  state 
of  stress  in  the  rock  mass  influences  the  stability  of  the  rock 
elements  or  discontinuities  after  the  excavation  and  the  magnitude 
of  the  forces  which  have  to  be  applied  to  the  excavated  boundaries 
to  simulate  the  creation  of  the  opening.  Therefore,  the  higher 
the  initial  stresses  are,  the  greater  the  stress  changes  and 
movements  are  in  the  rock  mass  surrounding  the  excavation. 

Two  types  of  measurements  using  overcoring  techniques  were  made 
at  the  Morrow  Point  Powerplant  to  determine  the  initial  state  of 
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TABLE  1  -  SUMMARY  OF  MATERIAL  PROPERTIES  USED  IN  ANALYSES,  MORROW  POINT  POWERPLANT 
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stress  (Bureau  of  Reclamation,  1965).  The  measurements  were 
made  before  the  excavation  in  the  6  ft  diameter  exploratory 
tunnel  at  distances  varying  from  4  ft  to  10  ft  behind  the 
excavated  surface  of  the  tunnel.  Test  results  showed  that 
the  in-situ  stresses  along  the  exploratory  tunnel  at  the 
powerplant  varied  from  1400  psi  to  600  psi  for  vertical  stress 
and  2700  psi  to  170  psi  for  horizontal  stress  indicating  that 
(1)  the  in-situ  stresses  varied  with  location,  (2)  the  vertical 
stress  may  be  greater  than  the  gravity  stress,  and  (3)  the  hori¬ 
zontal  stress  may  be  greater  than  the  vertical  stress.  Because 
of  the  limited  data  with  regard  to  the  magnitude  of  the  in-situ 
stresses  and  the  difficulties  associated  with  estimating  their 
values  in  the  rock  mass,  the  initial  state  of  stress  used  in  the 
analysis  was  obtained  by  a  gravity  turn-on  analysis  of  the  rock 
mass  without  the  opening.  These  analyses  were  conducted  to 
account  for  the  effects  of  the  steep  valley  wall  located  in  the 
vicinity  of  the  powerplant  excavation.  For  these  analyses,  the 
higher  values  of  the  Poisson's  ratio  (=  0.49)  was  assumed  for 
the  rock  to  simulate  the  likelihood  of  a  high  horizontal  stress. 

The  values  of  the  initial  stresses  which  have  to  be  applied  at 
excavated  boundary  to  simulate  the  excavation  were  obtained  by 
a  method  similar  to  that  suggested  by  Clough  and  Duncan  (1969). 
The  nodal  point  stresses  on  the  excavated  boundary  were  estimated 
from  the  stresses  in  the  surrounding  elements.  The  procedure  is 
described  in  Appendix  A. 
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Analysis  Procedures 

As  indicated  previously,  two  cases, using  high  and  low  normal 
and  shear  joint  stiffnesses  for  the  shear  zones  were  analyzed 
utilizing  the  computer  program  developed  in  Phase  1.  For  each 
case,  the  initial  stresses  for  each  element  were  first  obtained 
by  performing  a  gravity  turn-on  analysis.  The  nodal  stresses 
were  then  determined  along  the  excavated  surface  in  accordance 
with  the  procedure  described  in  Appendix  A  and  applied  at  the 
excavated  boundary.  Very  small  values  of  the  elastic  constants 
were  assigned  to  those  elements  in  the  opening  to  simulate  the 
cavity.  Material  properties  utilized  in  the  analysis  are  pre¬ 
sented  on  Table  1. 

Presentation  and  Discussion  of  Results 

The  results  of  the  analyses  are  presented  in  Figs.  45  and  46 
and  Table  2.  Fig.  45  illustrates  the  horizontal  displacements 
for  points  along  the  face  of  the  powerplant  chamber,  and  Fig.  46 
shows  the  vertical  displacements.  The  observed  movements  on 
the  a-  and  b-line  walls  are  also  shown  in  Figs.  45  and  46.  For 
the  case  with  high  joint  stiffnesses,  the  horizontal  inward 
movement  at  El.  6793  was  calculated  to  be  0.44  inch  on  the  a-line 
wall  and  0.41  inch  on  the  b-line  wall.  The  observed  inward  move¬ 
ments  at  the  corresponding  points  are  2.5  inches  on  the  a-line 
and  0.1  inch  on  the  b-line  wall.  When  the  low  joint  stiffnesses 
were  used,  the  calculated  inward  movements  increased  to  1.35  inches 
on  the  a-line  wall  and  0.62  inch  on  the  b-line  wall. 
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TABLE  2  -  OBSERVED  AND  CALCULATED  ROCK  MOVEMENTS  ON 
a-  AND  b- WALLS  IN  THE  POWERPLANT  CHAMBER 


Calculated*  (in.) 

Observed* 

(in.) 

Case  1 

With  High 

Joint  Stiffness 

a- line  wall 
at  El.  6793 

-2.5 

-0.44 

-1.35 

a-line  wall 
at  El.  6823 

-2.1 

-0.21 

-0.92 

b-line  wall 
at  El.  6793 

-0.1 

-0.41 

-0.62 

b-line  wall 
at  El.  6823 

+  0.2 

-0.18 

-0.37 

*  Only  horizontal  movements  are  indicated. 

»*-•»  indicates  movements  towards  centerline  of  chamber. 
"+"  indicates  movements  away  from  centerline  of  chamber. 
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The  results  indicate  that  the  analysis  provided  a  qualitative 
estimate  of  rock  movements  at  the  Morrow  Point  Powerplant 
excavation.  The  computed  deformations  are  of  the  same  order 
of  magnitude  as  those  observed.  The  difference  between  the 
computed  and  observed  deformation  is  significant  in  terms  of 
percentages  but  is  considered  satisfactory  in  terms  of  numer¬ 
ical  values  when  it  is  recognized  that  the  following  approxi¬ 
mations  were  made  in  the  analysis.  (1)  The  rock  wedge  and  the 
movements  associated  with  it  are  three-dimensional  in  nature, 
Fig.  42.  A  two-dimensional  plane  strain  analyses  would  tend 
to  underestimate  the  movements.  (2)  The  initial  state  of  stress 
was  estimated  in  the  analyses.  It  is  quite  possible  that  a 
higher  initial  state  of  stress,  especially  higher  horizontal 
stresses  might  exist  in  the  rock  mass  which  would  cause  higher 
induced  stresses  and  movements.  (3)  No  data  were  available  on 
the  properties  of  the  shear  zones  present  in  the  rock  mass.  It 
was  necessary  to  assume  their  values  in  the  analysis. 

It  should  also  be  recognized  that  the  powerplant  was  excavated 
in  stages,  whereas  the  analysis  assumed  that  the  excavation  was 
created  instantaneously.  For  an  ideal  elastic  material,  the 
assumption  would  not  result  in  any  error.  However,  in  those 
cases  where  the  loading  path  has  an  influence  on  the  results, 
the  effect  of  the  above  stated  assumption  can  be  significant. 
High  stress  concentrations  during  the  excavation  process  could 
cause  displacements  which  would  be  larger  than  those  computed 
on  the  assumption  of  an  instantaneous  creation  of  the  opening. 
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Summary  of  the  Achievement  for  Phase  2 

A  number  of  laboratory  studies  and  field  case  histories  were 
evaluated  to  determine  their  suitability  for  analysis  utilizing 
the  computer  program  developed  in  Phase  1.  A  laboratory  model 
study  of  an  opening  in  a  rock-like  material,  reported  by  Heuer 
and  Hendron  (1971) ,  and  the  excavation  for  the  Morrow  Point 
Powerplant  were  analyzed.  The  results  of  the  analyses  were 
compared  with  observed  behavior.  It  was  found  that  the  results 

of  the  analysis  could  be  utilized  to  predict  zones  of  potential 
failure . 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  puspose  of  this  study  was  to  evaluate  the  ability  of  the 
available  finite  element  techniques  for  the  solution  of  plane 
problems  to  predict  the  stresses,  strains,  and  deformations  in 
a  rock  mass  surrounding  an  excavation.  In  the  course  of  this 
investigation,  a  general  computer  program,  which  incorporates 
no  tension,  joint  perturbation  and  elasto-plastic  analysis,  was 
developed.  This  general  computer  program  was  employed  to  analyze 
the  Morrow  Point  Powerplant  excavation  and  a  laboratory  model 
study  of  an  opening  in  a  rock-like  material  conducted  by  Heuer 
and  Hendron  (1971).  The  results  of  the  analysis  were  compared 
with  actual  performance  . 

The  analysis  of  the  model  study  indicated  that  the  combined 
computer  program  could  identify  the  likely  failure  zones  that 
will  occur  around  an  opening  in  an  intact  rock.  However,  the 
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techniques  developed  are  not  adequate  to  simulate  accurately  the 
propagation  of  fractures  and  the  resulting  deformation  that  might 
occur. 

The  analyses  of  the  Morrow  Point  Underground  Powerplant  excavation 
indicated  that  the  developed  computer  program  could  be  utilized 
to  analyze  and  assist  in  predicting  the  performance  of  excavations 
in  a  rock  mass.  The  results  of  the  analyses  of  the  Morrow  Point 
Powerplant  excavation  provided  a  qualitative  estimate  of  the  move¬ 
ments  that  occurred  at  the  excavation.  The  accuracy  of  the  pre¬ 
diction  would  be  improved  if  better  information  was  available  on 
(1)  the  initial  state  of  stress  in  rock,  (2)  the  geologic  discon¬ 
tinuities  such  as  joints,  shear  zones,  foliations,  and  bedding 
planes,  and  (3)  in-situ  properties  of  the  rock. 

It  was  found  that  the  effort  necessary  to  idealize  problems  for 
analysis  was  far  greater  than  anticipated.  If  accurate  results 
are  desired,  great  care  has  to  be  taken  in  the  idealization. 

Since  real  problems,  in  general,  include  certain  factors  which 
cannot  be  modelled  accurately,  it  may  be  necessary  to  utilize 
more  than  one  idealization  to  obtain  bounds  on  the  possible 
behavior  of  the  excavation. 

Computational  techniques  which  have  the  capability  of  including 
realistic  idealizations  of  in-situ  rock  conditions  can  provide 
a  qualitative  estimate  of  the  behavior  of  a  rock  mass  in  the 
vicinity  of  an  excavation.  Considering  the  idealizations  used 
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and  the  accuracy  of  the  available  input  information,  the  ability 
of  the  computational  techniques  to  predict  quantitative  behavior 
is  considered  acceptable. 

Recommendations  for  Future  Research 

It  is  recommended  that  a  major  effort  should  be  directed  towards 
the  study  of  case  histories  for  the  purpose  of  establishing  the 
reliability  of  various  analytical  methods  to  predict  performance 
and  to  provide  a  basis  for  improvements.  Such  studies  should 
emphasize  the  requirements  for  idealization  and  field  input 
information.  Without  such  studies,  analytical  techniques  will 
not  gain  acceptance  in  the  design  profession. 

Ihis  study  has  provided  the  basis  for  recommending  that  the 

following  additional  features  be  incorporated  into  this  computer 
program: 

(i)  Theoretical  concepts  and  computational  techniques  for 
predicting  the  development  and  propagation  of  fractures 
in  a  rock  mass . 

(ii)  Computational  techniques  to  analyze  the  influence  of 
support  systems  and  construction  procedures. 

In  some  cases,  new  techniques  will  have  to  be  developed;  in 
others,  it  will  be  necessary  to  incorporate  work  already  done 
into  the  computer  program. 
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Operational  programs  for  conducting  a  linear  elastic  analysis 
of  three-dimensional  structures  are  available.  Efforts  should 
be  directed  towards  modifying  these  programs  to  include  the 
capability  of  modelling  and  analyzing  realistic  in-situ  rock 
conditions  e.g.,  joints,  no  tension,  etc.  Case  history  studies 
should  be  analyzed  to  determine  the  applicability  of  these  tech¬ 
niques  to  practical  problems. 
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APPENDIX  A 

EVALUATION  OF  NODAL  STRESSES  FROM  STRESSES 
IN  THE  SURROUNDING  ELEMENTS 

Dunlop,  Duncan  and  Seed  (1968),  Chang  and  Duncan  (1970),  and  Clough 
and  Duncan  (1969)  have  shown  that  excavation  may  be  simulated  in 
the  finite  element  method  by  applying  stresses  to  the  boundary  ex¬ 
posed  by  excavation.  This  technique  was  employed  in  the  present 
study.  In  the  finite  element  method,  stresses  are  generally  de¬ 
termined  only  at  the  element  centroids.  To  simulate  excavation,  it 
is  necessary  to  determine  the  stresses  on  the  surface  to  be  exposed 
by  excavation.  A  technique  similar  to  that  used  by  Clough  and  Duncan 
(1969)  was  employed  to  evaluate  nodal  stresses  on  the  excavated 
boundary  from  the  stresses  in  the  surrounding  elements. 

The  relationship,  between  the  known  stresses  at  the  element  centroids 
and  the  unknown  stresses  at  the  nodal  points  on  the  excavated  boun¬ 
dary  may  be  developed  using  an  interpolation  function  expressed  in 
the  following  form: 


a  »  a  +  a  x  +  a  y  +  a  xy  (A-l) 

1  2  3  4  ' 

in  which  a  is  the  nodal  stress  to  be  interpolated,  x  and  y  are  the 

’  •'  P 

coordinates  of  the  nodal  point  and  a  ,  a  ,  a  ,  and  a  are  inter1"" 

1  2  3  4 

polation  coefficients.  In  order  to  use  equation  (A-l)  to  determine 
the  stresses  at  the  nodes  of  an  element  to  be  excavated,  three 


WOODWARD-LUNDGREN  &  ASSOCIATES 


A-  3 


sets  of  the  interpolation  coefficients  are  calculated  (one  each 
f°r  crx,  0^,  and  n^)  using  the  known  stresses  in  the  element  to  be 
excavated  and  the  stresses  in  three  adjacent  elements.  For  a  given 
stress  component,  a,  the  unknown  interpolation  coefficients  are 
expressed  as: 


o(l)  -  a  +  a  x  +  a  y  +  a  x  y 

1  2  1  31  411 


o(2)  =  a  +  a  x  +ay  +  a  x  y 

1  2  2  3  2  4  2  2 

(A-  2) 

o(3)  -  *  a2X3  *  ^  *  Ws 


a(4)  =  a  +  ax  +  a  y  +  axy 

1  24  3  4  4  4  4 


in  which  a  (1)  is  the  stress  in  element  1,  a  (2)  is  the  stress  in 
element  2,  etc.  Equation  (A- 2)  may  be  expressed  in  matrix  form  as: 


{a}  =  [M]  {a}  (A- 3) 

in  which  (a)  is  the  stresses  for  elements  1,  2,  3  and  4,  [M]  is  the 
coordinate  matrix  for  elements  1,  2,  3  and  4,  and  {a}  is  the  unknown 
interpolation  coefficients,  {a}  may  be  expressed  by: 

{a}  =  [Mj " 1  io}a  CA-4) 

© 

The  interpolation  coefficients  may  then  be  used  to  solve  for  the 
stresses,  {a>n,  at  each  node  of  the  element  to  be  excavated. 
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(o}n  -  [N]  {a}  (A-5) 

in  which  [N]  is  the  coordinate  matrix  for  nodes  I,  J,  K,  and  L. 
Equations  (A-4)  and  (A-S)  may  be  combined,  as  follows: 

{a}  =  [N]  [M] " 1  {a>  (A-6) 

n  c 

Thus  the  values  of  stresses  at  the  nodes  of  an  element  to  be  ex¬ 
cavated  may  be  defined  in  terms  of  the  center  point  stresses  of 
that  element  and  three  adjacent  elements.  Evaluation  of  equation 
(A-6)  is  performed  in  Subroutine  NPSTRS, 
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A  Combined  Computer  Program  Using  Finite 
Element  Techniques  for  Elasto-plastic ,  Joint 
Perturbation,  and  No  Tension  Analysis  of 
Underground  Openings  in  Rock 
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A  COMBINED  COMPUTER  PROGRAM 
USING  FINITE  ELEMENT  TECHNIQUES  FOR  ELASTO- PLASTIC, 
JOINT  PERTURBATION,  AND  NO  TENSION  ANALYSIS 
OF  UNDERGROUND  OPENINGS  IN  ROCK 


Identification 

The  program  which  consists  of  a  main  program  and  13  subrou¬ 
tines  (STIFF,  MODIFY,  QUAD,  TRISTF ,  JTSTIF,  BANSOL,  STRESS,  REDO, 
JTSTR,  EPLAST,  INITST,  PRINST>  NPSTRS) ,  was  developed  by  Chin- 
Yung  Chang  using  a  computer  program  written  by  E.  Wilson  and 
modified  by  Goodman,  Taylor  and  Brekke  (1968). 

Purpose 

The  combined  computer  program  was  developed  on  the  basis  of 
three  concepts:  elasto-plastic ,  joint  perturbation  and  no  tension 
analysis  for  calculating  stresses  and  displacements  for  plane 
strain  conditions  in  a  rock  mass  surrounding  underground  openings. 

The  rock  mass  may  consist  of  joints,  faults,  bedding  planes  and 
other  geologic  discontinuities,  and  exhibit  elasto-plastic  and 
"no  tension"  behavior. 

Sequence  of  Operation 

(a)  The  main  program  handles  the  initial  input  and  monitors  the 
calling  of  the  subroutines  in  a  specified  order  as  shown  in 
Fig.  Al.  If  specified, for  the  last  iteration  of  the  last  incre¬ 
ment  in  an  analysis,  stresses  and  excess  stresses  to  be  re¬ 
distributed  for  two-dimensional  elements,  normal  and  tangential 
stresses  for  joint  elements,  nodal  point  displacements  and  yield 
functions  for  two-dimensional  elements  are  punched  onto  cards 
to  be  used  for  restarting  computation. 


(b)  Subroutine  STIFF  assembles  the  general  stiffness  matrix  for 

the  entire  structure,  adds  in  concentrated  loads  at  the  nodal 
points,  adds  in  loads  due  to  boundary  pressures  and  modifies 
the  stiffness  matrix  for  the  boundary  conditions. 


(c)  In  Subroutine  QUAD  the  constitutive  equations  are  formulated. 

(d)  Subroutine  TRISTF  forms  the  stiffness  matrix  for  triangular 
subelements  and  if  specified*,  element  loadings  due  to  gravity 
are  generated. 


(e)  Subroutine  JTSTIF  forms  the  stiffness  matrix  for  each  joint 
element . 
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FIG.B1  -  SIMPLIFIED  FLOW  DIAGRAM  SHOWING  SEQUENCE  OF  OPERATION  OF  ALL  SUBROUTINES 
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(f)  Subroutine  MODIFY  modifies  the  stiffness  matrix  for 
the  boundary  conditions. 


(a)  Subroutine  REDO  calculates  equilibrating  nodal  point 
forces  due  to  gravity  if  specified,  and  for  excess 
stresses  if  the  elements  are  subjected  to  tensile 
stresses  greater  than  that  allowable. 


fh)  Subroutine  BANSOL  solves  the  simultaneous  equations 
representing  the  structural  stiffness  matrix  and 
the  structural  load  vector  for  nodal  point  dis¬ 
placements  . 

m  Subroutine  STRESS  calculates  incremental  stresses  and 
strains,  cumulates  stresses,  and  prints  stresses  and 
strains  for  two-dimensional  elements. 


(j)  Subroutine  JTSTRS  calculates  and  prints  normal  and 
tangential  stresses,  normal  and  tangential  dis¬ 
placements  (cumulative  and  incremental)  and  excess 
normal  and  tangential  stresses  to  be  redistributed 
by  comparing  stress  with  strength  for  joint  elements. 
The  equilibrating  nodal  point  forces  are  also  com- 
puted  from  the  excess  stresses  and  stored  for  the 
next  iteration. 

(k)  Subroutine  EPLAST  calculates  yield  functions  and 
elasto-plastic  stress-strain  relation  for  those 
two-dimensional  elements  in  yield.  The  excess 
stresses  to  be  redistributed  are  computed  as  a 
difference  between  changes  in  stress  calculated 
from  the  elastic  stress -strain  relation  and  those 
calculated  from  the  elasto-plastic  stress -strain 
relation. 

(l)  Subroutine  INITST  generates  initial  stresses  under 
free-field  conditions. 

(m)  Subroutine  PRINST  calculates  magnitudes  and  directions 
of  the  principal  stresses  and  strains. 

(n)  Subroutine  NPSTRS  computes  stresses  at  nodes  on  the  ex¬ 
cavated  boundary  from  stresses  in  surrounding  elements. 
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Output 

The  data  describing  the  finite  element  configuration, 
the  material  properties  and  pressures  applied  to  the  exca¬ 
vated  face  to  simulate  excavation  for  the  opening  are  printed 
after  being  read.  Nodal  point  displacements  (incremental  and 
cumulative),  stresses,  strains  and  yield  functions  for  two- 
dimensional  elements;  normal  and  tangential  stresses,  normal 
and  tangential  displacements  (incremental  and  cumulative)  for 
joint  elements,  are  printed  after  each  increment  or  iteration. 
If  specified,  for  the  last  iteration  of  the  last  increment 
in  an  analysis,  stresses,  strains  and  excess  stresses  to  be 
redistributed  for  two-dimensional  elements,  normal  and  tan¬ 
gential  stresses  for  joint  elements,  nodal  point  displacements 
and  yield  functions  for  two-dimensional  elements  are  punched 
onto  cards  to  be  used  for  restarting  computation. 

Input  Data  Procedure 

1st  CARD  TYPE:  FORMAT  (8A10)  (One  Card) 

Cols  2-80  Identifying  information  to  be  printed  with 
results . 


2nd  CARD  TYPE: 

FORMAT 

(415,  2F10.2,  215,  3F10.5) 

(One  Card) 

Cols.  1-5 

NUMNP  - 

Number  of  nodal  points 
(maximum  900) 

6-10 

NUMEL  - 

Number  of  elements 
(maximum  800) 

11-15 

NUMMAT  - 

Number  of  different  materials 
(maximum  12) 

16-20 

NUMPC  - 

Number  of  pressure  cards 

21-30 

ACELX  - 

Acceleration  in  X-direction 

31-40 

ACELY  - 

Acceleration  in  Y-direction 

41-45 

NP 

Number  of  approximations  (increments) 

46-50 

NRES*  - 

=  -1,  Residual  stresses  generated 
from  which  residual  load  is 
calculated. 

*IP  NREAD  =  1,  NRES  should  not  be  greater  than  zero. 

';f  NRES  =  -1  or  1,  gravity  turn-on  analysis  is  performed. 
If  NRES  =  -1„  stresses  at  nodes  on  the  excavated  boundary 
are  to  be  calculated  and  punched  out  onto  cards,  i.e. 

2 1th1  and  22nd  card  types  are  needed. 
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=  0,  Residual  stresses  generated 
but  residual  load  is  zero. 

=  l,  Residual  stresses  read  as  in¬ 
put  from  which  residual  load 
is  generated. 

=  2,  Residual  stresses  read  as  in¬ 
put,  but  residual  load  is  zero. 

Cols.  51-60  FRAC  -  Percentage  of  maximum  tensile  stress 

considered  as  cracked  zone,  read  in 

zero . 

61-70  REFPR  -  Vertical  stress  at  the  reference 

point. 

71-80  DEPTH  -  Y  -  ordinate  at  the  reference  point. 

3rd  CARD  TYPE:  FORMAT  (1615)  (One  Card) 

Cols.  1-5  NPRSNT  -  Present  loading  increment  number. 

6-10  NREAD  -  =  0,  no  data  from  previous  computa¬ 

tion  will  be  read  as  input. 

=  1,  data  from  last  increment  are 
read  as  input. 

11-15  NPUNCH  -  =■  0,  data  will  not  be  punched  out 

at  the  last  iteration. 

=  1,  data  will  be  punched  out  at  the 
last  iteration  of  the  last  increment. 

16-20  NSTSRT  -  ■  1,  stresses  in  R-0  directions  will 

be  printed.  For  this  case,  the 
center  of  the  circular  opening 
has  to  be  located  at  the  origin 
of  the  coordinate  system. 

4th  CARD  TYPE:  FORMAT  (15,  F10.5) 

Cols.  1-5  ITN(N)  -  Number  of  iterations  for  N 

increment. 

6-15  PRATIO(N)  -  Percentage  of  total  pressure  ap¬ 
plied  for  Nth  increment. 

Repeat  for  each  loading  increment. 

5th  CARD  TYPE:  FORMAT  (1615)  (One  Card) 

Cols  1-5  MJ0INT  -  Total  number  of  material  types  for 

joints  (maximum  12) 

6-10  MTENS  -  Total  number  of  material  types  that 

can  sustain  tension. 
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6th  CARD  TYPE:  FORMAT  (1615)  (Omit  this  card  if  MJOINT=0 

on  5th  card  type) 

Cols.  1-5  MJNT(I)  -  Material  type  number  for  joint  elements. 

6-10  Same 
11-15  - 

7th  CARD  TYPE:  FORMAT  (1615)  (Omit  this  card  if  MTENS=0 

on  5th  card  type) 

Cols.  1-5  MNTEN(I)  -  Material  type  number  which  can  sustain 

tension. 

6-10  Same 
11-15  - 

8th  CARD  TYPE:  FORMAT  (15,  2F10.5)  (One  Card) 

Cols.  1-5  MTYPE  -  Material  type  number. 

6-15  RO (MTYPE)-  Mass  density  of  this  material  type. 

16-20  AKO (MTYPE)  -  Ratio  of  horizontal  to  vertical 

stress. 

9th  CARD  TYPE:  FORMAT  (8F10.5) 

Cols.  1-10  E (1 ,  MTYPE)-  Tensile  strength  for  normal  materials 

or  normal  stiffness  for  joint  materials 

11-20  E(2,  MTYPE)-  Modulus  in  compression  for  normal  mtls. 

or  shear  stiffness  for  joint  materials 

21-30  E (3 ,  MTYPE)-  Poisson's  ratio  for  normal  materials 

or  cohesion  for  joint  materials 

31-40  E (4 ,  MTYPE)-  Modulus  in  tension  for  normal  materials 

angle  of  friction  for  joint  mtls. (degree 

41-50  E(5,  MTYPE)-  Cohesion  for  normal  materials  or 

maximum  allowable  closure  (input  as 
negative)  for  joint  materials 

51-60  E (6 ,  MTYPE)-  Angle  of  friction  for  normal  mtls. 

(degrees) 

Repeat  8th  and  9th  card  types  for  all  material  types. 

10th  CARD  TYPE:  FORMAT  (15,  F5.0,  4F10.0)  (One  card  for  each 

nodal  point) 

Cols.  1-5  N  -  Nodal  point  number 

6.-10  CODE  (N)-  Number  which  indicates  if  displacements 

or  forces  are  to  be  specified 

=0  UR  is  the  specified  X-load  and 
UZ  is  the  specified  Y-load 
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=1  UR  is  the  specified  X-displacement 
and  UZ  is  the  specified  Y-load 

=2  UR  is  the  specified  X-load  and 

UZ  is  the  specified  Y-displacement 

=3  UR  is  the  specified  X-displacement 
and  UZ  is  the  specified  Y-displace¬ 
ment 


Cols.  11-20 

R(N) 

X-ordinate 

21-30 

Z(N) 

Y-ordinate 

31-40 

UR(N) 

X-load  or  displacement 

41-50 

UZ  (N) 

Y-load  or  displacement 

Nodal  points  must  be  numbered  in  sequence.  If  nodal  point  numbers 
are  omitted,  those  omitted  are  generated  automatically  at  equal 
spacings,  between  those  specified  and  CODE(N)  is  assigned  zero. 

The  first  and  last  must  be  specified. 


11th  CARD  TYPE:  FORMAT  (1615)  (One  card  for  each  element) 


1-5 

M 

6-10 

IX(M,1)- 

11-15 

IX (M, 2)  - 

16-20 

IX (M, 3) - 

21-25 

IX(M,4)  - 

26-30 

IX(M,5) - 

Element  number 
Nodal  point  I 
Nodal  point  J 
Nodal  point  K 
Nodal  point  L 
Material  number 


The  nodal  point  numbers  must  be  numbered  consecutively  proceeding 
counterclockwise  around  the  elements.  The  nodal  point  numbers 
for  any  element  must  not  differ  by  more  than  39.  If  element  num¬ 
bers  are  omitted,  those  missing  will  be  generated  by  incrementing 
the  element  number  and  each  nodal  point  number  (I,  J,  K  ana  w 
by  one,  and  assigning  the  same  material  number  as  the  last  element 
specified.  The  first  and  last  elements  must  be  specified. 

Triangular  elements  are  also  permissable,  and  are  identified  by 
repeating  the  last  nodal  point  number  (i,e.,  I,  J,  K,  K) . 
ioint  elements,  nodal  points  must  be  numbered  I,  J,  K,  L  counter 

clockwise  proceeding  along  length  of  j®*-?1  from  1  **  f  r 

length  from  K  to  L.  Nodal  points  I  and  L  (J  and  K)  have  different 

numbers  but  identical  coordinates. 


One -dimensional  elements  are  also  permissable  and  are  identified 
by  the  node  number  sequence  (I,  J,  J,  I) • 
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12th  CARD  TYPE;  FORMAT  (1615) 

Cols.  1-5  I  -  Nodal  point  number  I  along  the  boundary  IJ 

where  the  boundary  pressure  is  applied. 

6-10  J  -  Nodal  point  number  J  along  the  boundary  IJ. 

11-15)  Same  as  above;  two  nodal  point  numbers  for 

- h  each  boundary  where  the  boundary  pressure 

16~20J  is  applied. 

21-25! 

26  -30; 

As  shown  in  Figure  B2,  nodal  points  I  and  J  must  be  ordered 
in  counterclockwise  order  about  centroid  of  element  on  which 
the  pressure  is  applied. 

15th  CARD  TYPE:  FORMAT  (1615) 

Cols.  1-5  NPCAV  -  Number  of  nodal  points  along  the  boundary 

where  the  boundary  pressures  are  applied. 

14th  CARD  TYPE:  FORMAT  (15,  3F10.0) 

Cols.  1-5  NPCA(M)  -  Nodal  point  number  where  the  boundary 

pressure  is  applied. 

6-15  PSCA(M, 1)  -  X  normal  stress  at  nodal  point  NPCA(M) 

16-25  PSCA(M, 2)  -  Y  normal  stress  at  nodal  point  NPCA(M) 

26-35  PSCA (M, 3)  -  XY  shear  stress  at  nodal  point  NPCA(M) 


Stresses  (ox,  oyS 


and  t  )  are  shown  positive  in  Figure  B2 
xy 


As  shown  in  Figure  B2,  stresses  (ax,  ay>  xxy)  at  the  nodal 
point  are  converted  to  normal  and  shear  stress  on  the  boun¬ 
dary.  Then  the  nodal  point  forces  are  calculated  from  the 
normal  and  shear  stress  along  the  boundary  assuming  linear 
stress  distribution  along  the  boundary.  The  normal  and  shear 
stress  along  the  boundary  are  shown  positive  in  Figure  B2. 


12th,  13th  and  14th  card  types  are  neglected  if  NUMPC  =  0 
on  the  2nd  card  type. 


V 


15  CARD  TYPE;  FORMAT  (15,  3E15.5)  (One  card  for  each  element) 

Cols.  1-5  N  -  Element  number 

6-20  STRS (N, 1)  -  Initial  stress  ax 
21-35  STRS (N , 2)  -  Initial  stress 
36-50  STRS (N , 3)  -  Initial  stress  x 

xy 

This  card  type  is  neglected  if  NRES  <_  0. 
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FIG.  B2  -  SIGN  CONVENTION  FOR  BOUNDARY  PRESSURE 
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16th  CARD  TYPE:  FORMAT  (15) 

Cols*  1-5  NJT  -  Total  number  of  joint  elements. 

17th  CARD  TYPE:  FORMAT  (2F20.5,  15) 

Cols.  1-20  FN(I)  -  Normal  stress  across  joint  element  NEJT(I) 

21-40  FT(I)  -  Tangential  stress  across  joint  element 
NEJT(I) 

41-45  NEJT(I)-  Element  number  for  joint  element 

Repeat  17th  card  type  for  all  joint  elements 

16th  and  17th  card  types  are  neglected  if  NRES  £  0  on  the 
2nd  card  type. 


18th  CARD  TYPE:  Binary  data  cards  for  all  elements 


STRS(N,1)  -  ax 
STRS(N,2)  -  ov 
STRS(N,5)  -  txy 


SEP(N,1)  -  Excess  stress  ox* 

SEP (N,  2)  -  Excess  stress  ay’ 

SEP (N,  3)  -  Excess  stress  Tx'y 

MTAG(N)  -  An  index  .indicating  if  the  element  has  failed  in 
tension. 

If  MTAG(N)  »  1,  the  element  has  not  failed  in 
tension. 

MTAG(N)  =  2,  the  element  has  failed  in  tension 
in  one  direction  only 

MTAG(N)  =  3,  the  element  has  failed  in  tension 
in  two  orthogonal  directions. 


19th  CARD  TYPE:  Binary  data  cards  for  all  joint  elements, 

neglected  if  MJOINT  =  0. 


NJT  -  Total  number  of  joint  elements 
FN(N)  -  Normal  stress  across  joint  element 
FT(N)  -  Tangential  stress  across  joint  element 

18th  and  19th  card  types  are  ngelected  if  NREAD  *  0  on  the  3rd 
card  type,  i.e.,  these  cards  are  need  to  restart  the  computation. 

20th  CARD  TYPE:  Binary  data  cards  for  all  nodal  point  dis¬ 
placements,  yield  function,  axial  stress,  and 
all  strain  components  for  all  elements. 

DISP(N,1)  -  x  -  displacement  for  node  N 
DISP(N,2)  -  y  -  displacement  for  node  N 

FY(N)  -  Yield  function  for  element  N 
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SIGZ(N)  -  axial  stress  for  element  N 


STRN (N, 1)  -  strain  component,  sx,  for  element  N 
STRN (N, 2)  -  strain  component,  ey,  for  element  N 
STRN (N, 3)  -  strain  component,  exy,  for  element  N 

20th  card  type  is  neglected  if  NREAD  =  0  on  the  3rd  card 
type,  i.e.,  these  cards  are  needed  to  restart  the  computation. 

21th  CARD  TYPE:  FORMAT  (15) 

Cols.  1-5  NPST  -  Number  of  nodal  points  at  which  nodal  stresses 

are  to  be  computed. 

22nd  CARD  TYPE;  FORMAT  (1615) 


Cols.  1-5 


6-10 1 
11-15  I 
16-20 
21-25 
26-30 
31-351 
36-40  [ 
41-45 
46-50 


NS(I)  -  Nodal  point  number  at  which  nodal  stresses 
are  to  be  computed. 

NSEL(I,1)  -  First  interpolation  element  number 
NSEL(I,2)  -  Second  interpolation  element  number 
NSEL(I,3)  -  Third  interpolation  element  number 
NSEL(I,4)  -  Fourth  interpolation  element  number 


same  as  above 


Repeat  for  all  nodal  points  at  which  nodal  stresses  are  to 
be  computed.. 

The  mid-points  of  no  three  of  the  four  interpolation  ele¬ 
ments  may  lie  on  a  horizontal  or  vertical  line.  These  ele¬ 
ments  should  be  read  in  a  criss-crossed  fashion  as  shown  in 
Fig.  B3.  The  centroids  of  the  first  and  second  elements 
should  not  lie  on  a  vertical  line. 

21th  and  22nd  card  types  are  required  if  NRES  *  -1  on  the  2nd 
card  type. 
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FIG.  B3  -  SEQUENCE  FOR  READING  IN 
INTERPOLATION  ELEMENTS 


is?. 
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CAPABILITIES  OF  OPERATIONAL  THREE-DIMENSIONAL 
FINITE  ELEMENT  PROGRAMS  FOR  THE  ANALYSIS 
OF  PROBLEMS  IN  ROCK  MECHANICS 


Introduction 


With  very  few  exceptions,  all  the  work  done  in  the  application  of 
finite  element  methods  to  problems  in  rock  mechanics  has  been  con¬ 
cerned  with  axisymmetric  and  plane  problems.  In  general,  an  actual 
field  situation  is  three-dimensional  and,  in  some  cases,  an  axi¬ 
symmetric  or  plane  formulation  of  the  problem  may  not  be  satisfac¬ 
tory.  Therefore,  the  development  of  finite  element  programs  for 
the  three-dimensional  analysis  of  boundary  value  problems  is  of 
interest  in  design  of  excavations  in  rock. 

The  purpose  of  this  brief  discussion  is  not  to  elaborate  on  the 
theory  or  the  details  of  the  development  of  three-dimensional 
finite  element  programs,  but  to  evaluate  what  is  available  in 
terms  of  operational  programs  in  their  ability  to  assist  in  the 
design  of  excavations  in  rock. 

Capabilities  of  Operational  Programs 

When  the  original  proposal  was  written,  there  was  very  little  avail¬ 
able  in  documentation  on  operational  programs.  Since  that  time, 
some  information  has  become  available  on  existing  programs.  The 
most  significant  piece  of  information  has  been  a  recent  report  by 
Professor  E.  L.  Wilson*,  which  provides  detailed  information  on  an 
operational  program  for  the  analysis  of  three-dimensional  solid 
structures.  Since  this  is  the  most  recent  reference  available,  it 
is  utilized  as  a  basis  for  considering  the  available  capabilities 
for  analyzing  excavations  in  rock. 

It  is  important  to  recognize  that,  as  in  the  historical  development 
of  plane  finite  element  programs  for  use  in  practical  problems  in 
rock  mechanics,  the  first  step  is  the  development  of  the  program 
to  analyze  linear  elastic  problems.  From  a  fundamental  standpoint, 
this  is  the  most  difficult  step.  The  attainment  of  this  is  there¬ 
fore  of  considerable  significance.  It  is  now  a  question  of  adapt¬ 
ing  this  basic  program  as  was  done  for  the  plane  case  to  include 
the  various  features  that  are  significant  in  rock  mechanics  prob¬ 
lems.  This  is  no  trivial  task,  but  the  most  important  step  has 
been  taken. 

For  examining  the  capabilities  of  the  available  program,  it  is 
most  convenient  to  consider  the  factors  that  have  to  be  modelled 
in  evaluating  the  stability  of  openings  in  rock  and  examine  the 


*"SOLID  SAP,  A  Static  Analysis  Program  for  Three-Dimensional  Solid 
Structures,  UC  SESM  71-19"  by  Edward  L.  Wilson.  Report  to  Denver 
Mining  Research  Center,  U.S.  Department  of  Interior,  Structural 
Engineering  Laboratory,  University  of  California,  Berkeley, 
California.  September  1971,  There  was  an- earlier  version  which 
was  not  as  suitable  for  rock  mechanics  problems  as  the  above  report. 
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ability  o£  the  available  finite  element  program  to  model  these 
factors . 

Geometry  -  In  general,  the  geometry  of  the  problem  is  three- 
dimensional  .  The  most  significant  capability  of  the  program  is 
its  ability  to  model  this  aspect  of  the  problem. 

Material  Properties  -  At  the  present  time,  the  computer  program 
only  considers  linear  elastic  materials.  The  three-dimensional 
elements  also  have  the  further  limitation  of  isotropy.  There 
are  no  provisions  for  conducting  a  "no  tension"  or  "elasto- 
plastic"  analysis  or  considering  any  other  nonlinear  material 
behavior. 

Geologic  Conditions  -  The  two  major  considerations  are  non- 
homogeneity  as  evidenced  by  layered  deposits  and  the  existence 
of  joints  and  other  discontinuities.  The  non-homogeneity  of 
the  rock  can  be  readily  modelled.  However,  there  are  no  pro- 
visions  for  considering  joints  and  discontinuities. 

An  earlier  published  paper  by  Mahtab  and  Goodman*  discusses  the 
three-dimensional  finite  element  analysis  of  jointed  rock  slopes. 
This  paper  is  based  on  a  doctoral  dissertation  by  Mahtab.**  The 
use  of  the  program  was  demonstrated  by  solving  a  special  case  of 
a  symmetric  rock  wedge  on  two  joint  planes  which  had  linear 
stress-deformation  characteristics.  The  general  capability  of  the 
program  is  not  clear.  It  does  appear  that  the  program  developed 
by  Mahtab  does  not  have  the  computational  refinements,  the  number 
of  element  types,  the  superior  element  characteristics  and  docu¬ 
mentation,  but  it  does  have  the  capability  of  modelling  joints. 

Construction  Sequence  -  Since  the  capabilities  of  the  program 
are  confined  to  linear  elastic  systems,  the  construction  sequence 
has  no  significance;  hence,  there  are  no  provisions  for  con¬ 
sidering  this  factor. 

Support  Schemes  -  The  program  has  the  capability  for  analyzing 
various  support  schemes.  The  different  element  types  e.g., 
plate,  shell,  truss,  beam,  etc.,  can  be  used  to  model  a  variety 
of  support  schemes.  However,  as  far  as  we  are  aware,  this  fea¬ 
ture  of  the  program  has  not  been  exploited. 

Except  for  the  example  solved  by  Professor  Wilson  and  presented  in 
the  above-mentioned  report,  there  does  not  appear  to  be  any  sig¬ 
nificant  information  on  the  computer  time  and  time  required  to 


*  Mahtab,  M.  A.  and  Goodman,  R.  E.,  "Three-Dimensional  Finite 
Element  Analysis  of  Jointed  Rock  Slopes,"  Procs.  2nd  Congress 
of  the  International  Society  of  Rock  Mechanics.  Vol  3,  No.  7- 
12,  Belgrade,  Yugoslavia,  1970. 

**  Mahtab,  M.  A.,  "Three  Dimensional  Finite  Element  Analysis  of 
Jointed  Rock  Slopes,"  Ph.D.  Thesis,  University  of  California, 
Berkeley,  1970. 
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prepare  the  data.  The  problem  solved  was  a  quarter  of  a  pillar 
and  room  in  a  mine  of  a  homogeneous  medium.  It  was  a  relatively 
simple  problem.  It  took  two  man  days  to  prepare  the  input  data. 
The  actual  computer  cost  on  CDC  6400  was  §50.00  at  commercial 
rates. 

Idealizations  would  take  considerably  more  time  and  the  computa¬ 
tional  costs  would  be  far  greater  when  there  are  joints  and  dif¬ 
ferent  materials  located  in  geometrically  irregular  patterns. 
Furthermore,  if  the  materials  and  joints  have  nonlinear  char¬ 
acteristics  and  the  program  is  modified  to  model  this  effect, 
the  computational  time  would  be  further  increased.  Once  itera¬ 
tion  and  convergence  questions  become  included,  the  time  and 
cost  in  solving  practical  problems  increases.  It  is  necessary 
to  evaluate  if  the  additional  cost  and  effort  in  idealizing  the 
problem  and  conducting  the  analysis  is  justified  in  terms  of  a 
better  prediction  of  performance.  This  can  only  be  done  on  the 
basis  of  comparing  predictions  and  observed  performance  for 
actual  field  problems. 

Final  Remarks 


It  would  appear  that  two  significant  steps  have  been  taken  (i)  a 
general  linear  elastic  program  has  been  developed  and  (ii)  the 
formulation  of  three-dimensional  joint  elements  has  been  completed. 
The  next  step  is  to  extend  the  capabilities  of  the  linear  elastic 
program  to  include  factors  that  are  necessary  to  realistically 
model  actual  field  problems.  Initially,  the  joint  elements  that 
have  been  developed  could  be  incorporated  into  the  program.  As  a 
further  refinement,  nonlinear  joint  behavior  which  is  controlled 
by  the  friction  and  cohesion  along  the  joint  should  be  included. 

Once  this  is  done,  problems  of  practical  interest  could  be  analyzed. 

It  is  also  necessary  to  gain  experience  in  the  modelling,  analyses 
and  interpretation  of  the  results  on  the  basis  of  case  history 
studies.  Some  of  this  information  should  be  obtained  in  an  actual 
design  environment  where  the  constraints  of  time,  money  and  avail¬ 
able  input  information  are  always  present. 


WOODWARD-LUNDGREN  &  ASSOCIATES 


APPENDIX  D 


TANGENTIAL  STRESS  DISTRIBUTIONS  AND 
DISPLACEMENTS  IN  ELASTO-PLASTIC  ANALYSIS 
OF  A  CIRCULAR  OPENING 


WOODWARD-LUNDGREN  &  ASSOCIATES 


TANGENTIAL  STRESS  DISTRIBUTIONS  AND 
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In  the  semi-annual  technical  report,  a  comparison  of  the  dis¬ 
tribution  of  tangential  stresses  obtained  from  an  elasto-plastic 
analysis  of  the  problem,  presented  in  Figure  Dl,  conducted  by 
Reyes  (1966),  Baker  et  al.  (1969),  and  the  study  under  Contract 
H0210046  was  presented.  The  results  of  this  comparison  are 
shown  in  Figures  D2  and  D3„  It  can  be  observed  that  there  are 
significant  differences  in  the  tangential  stress  distribution 
and  displacements.  The  differences  could  have  occurred  from 
three  sources:  (i)  the  different  loading  paths  used,  (ii) 
variations  in  the  method  of  computing  the  axial  stress,  and 
(iii)  the  degree  of  convergence  obtained. 

The  problem  was  reanalyzed  following  a  loading  path  consistent 
with  that  used  by  Baker,  et  al.  (1969),  and  similar  to  that 
used  by  Reyes  (1966) ,  by  applying  the  pressures  on  the  outer 
boundary.  The  results  of  the  analysis  are  shown  in  the  body 
of  the  report,  Figures  16  and  17.  Th^  stress  distribution 
and  the  peak  tangential  stress  obtained  are  in  satisfactory 
agreement  with  the  results  of  Baker  et  al.  and  Reyes,  differ¬ 
ences  being  of  the  order  of  +151,  and  the  displacements  are 
almost  identical.  The  results  presented  in  the  semi-annual 
technical  report  were  obtained  by  assigning  initial  stresses 
to  elements  surrounding  the  opening  and  applying  the  boundary 
pressure  to  the  excavated  boundary  to  simulate  the  creation 
of  the  opening.  The  difference  in  the  results,  when  compared 
with  those  obtained  in  the  reanalyses,  indicates  the  effect  of 
the  loading  or  unloading  path  on  the  stress  distribution  in  an 
elasto-plastic  material.  This  points  out  the  importance  of 
modelling  the  loading  sequence. 

Some  difference  in  the  tangential  stress  distribution  may  also 
be  attributed  to  the  degree  of  convergence  obtained  by  each 
method  of  analysis  and  the  different  methods  used  in  calculating 
the  axial  stress.  This  has  been  discussed  in  the  report. 
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FIG.  D2  -  VERTICAL  AND  HORIZONTAL  STRESSES  ALONG  HORIZONTAL 
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FIG.  D3  -  DEFORMATION  ALONG  CAVITY  FACE  OF  A  CIRCULAR  OPENING  AS 
COMPUTED  BY  ELASTIC  AND  ELASTIC-PLASTIC  ANALYSIS 


